text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_events_data
using AWS.Compat
using AWS.UUIDs

"""
    batch_acknowledge_alarm(acknowledge_action_requests)
    batch_acknowledge_alarm(acknowledge_action_requests, params::Dict{String,<:Any})

Acknowledges one or more alarms. The alarms change to the ACKNOWLEDGED state after you
acknowledge them.

# Arguments
- `acknowledge_action_requests`: The list of acknowledge action requests. You can specify
  up to 10 requests per operation.

"""
batch_acknowledge_alarm(acknowledgeActionRequests; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/acknowledge", Dict{String, Any}("acknowledgeActionRequests"=>acknowledgeActionRequests); aws_config=aws_config)
batch_acknowledge_alarm(acknowledgeActionRequests, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/acknowledge", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("acknowledgeActionRequests"=>acknowledgeActionRequests), params)); aws_config=aws_config)

"""
    batch_disable_alarm(disable_action_requests)
    batch_disable_alarm(disable_action_requests, params::Dict{String,<:Any})

Disables one or more alarms. The alarms change to the DISABLED state after you disable them.

# Arguments
- `disable_action_requests`: The list of disable action requests. You can specify up to 10
  requests per operation.

"""
batch_disable_alarm(disableActionRequests; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/disable", Dict{String, Any}("disableActionRequests"=>disableActionRequests); aws_config=aws_config)
batch_disable_alarm(disableActionRequests, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/disable", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("disableActionRequests"=>disableActionRequests), params)); aws_config=aws_config)

"""
    batch_enable_alarm(enable_action_requests)
    batch_enable_alarm(enable_action_requests, params::Dict{String,<:Any})

Enables one or more alarms. The alarms change to the NORMAL state after you enable them.

# Arguments
- `enable_action_requests`: The list of enable action requests. You can specify up to 10
  requests per operation.

"""
batch_enable_alarm(enableActionRequests; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/enable", Dict{String, Any}("enableActionRequests"=>enableActionRequests); aws_config=aws_config)
batch_enable_alarm(enableActionRequests, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/enable", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("enableActionRequests"=>enableActionRequests), params)); aws_config=aws_config)

"""
    batch_put_message(messages)
    batch_put_message(messages, params::Dict{String,<:Any})

Sends a set of messages to the AWS IoT Events system. Each message payload is transformed
into the input you specify (\"inputName\") and ingested into any detectors that monitor
that input. If multiple messages are sent, the order in which the messages are processed
isn't guaranteed. To guarantee ordering, you must send messages one at a time and wait for
a successful response.

# Arguments
- `messages`: The list of messages to send. Each message has the following format: '{
  \"messageId\": \"string\", \"inputName\": \"string\", \"payload\": \"string\"}'

"""
batch_put_message(messages; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/inputs/messages", Dict{String, Any}("messages"=>messages); aws_config=aws_config)
batch_put_message(messages, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/inputs/messages", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("messages"=>messages), params)); aws_config=aws_config)

"""
    batch_reset_alarm(reset_action_requests)
    batch_reset_alarm(reset_action_requests, params::Dict{String,<:Any})

Resets one or more alarms. The alarms return to the NORMAL state after you reset them.

# Arguments
- `reset_action_requests`: The list of reset action requests. You can specify up to 10
  requests per operation.

"""
batch_reset_alarm(resetActionRequests; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/reset", Dict{String, Any}("resetActionRequests"=>resetActionRequests); aws_config=aws_config)
batch_reset_alarm(resetActionRequests, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/reset", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resetActionRequests"=>resetActionRequests), params)); aws_config=aws_config)

"""
    batch_snooze_alarm(snooze_action_requests)
    batch_snooze_alarm(snooze_action_requests, params::Dict{String,<:Any})

Changes one or more alarms to the snooze mode. The alarms change to the SNOOZE_DISABLED
state after you set them to the snooze mode.

# Arguments
- `snooze_action_requests`: The list of snooze action requests. You can specify up to 10
  requests per operation.

"""
batch_snooze_alarm(snoozeActionRequests; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/snooze", Dict{String, Any}("snoozeActionRequests"=>snoozeActionRequests); aws_config=aws_config)
batch_snooze_alarm(snoozeActionRequests, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/alarms/snooze", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("snoozeActionRequests"=>snoozeActionRequests), params)); aws_config=aws_config)

"""
    batch_update_detector(detectors)
    batch_update_detector(detectors, params::Dict{String,<:Any})

Updates the state, variable values, and timer settings of one or more detectors (instances)
of a specified detector model.

# Arguments
- `detectors`: The list of detectors (instances) to update, along with the values to update.

"""
batch_update_detector(detectors; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/detectors", Dict{String, Any}("detectors"=>detectors); aws_config=aws_config)
batch_update_detector(detectors, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/detectors", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("detectors"=>detectors), params)); aws_config=aws_config)

"""
    describe_alarm(alarm_model_name)
    describe_alarm(alarm_model_name, params::Dict{String,<:Any})

Retrieves information about an alarm.

# Arguments
- `alarm_model_name`: The name of the alarm model.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"keyValue"`: The value of the key used as a filter to select only the alarms associated
  with the key.
"""
describe_alarm(alarmModelName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/alarms/$(alarmModelName)/keyValues/"; aws_config=aws_config)
describe_alarm(alarmModelName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/alarms/$(alarmModelName)/keyValues/", params; aws_config=aws_config)

"""
    describe_detector(detector_model_name)
    describe_detector(detector_model_name, params::Dict{String,<:Any})

Returns information about the specified detector (instance).

# Arguments
- `detector_model_name`: The name of the detector model whose detectors (instances) you
  want information about.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"keyValue"`: A filter used to limit results to detectors (instances) created because of
  the given key ID.
"""
describe_detector(detectorModelName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)/keyValues/"; aws_config=aws_config)
describe_detector(detectorModelName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)/keyValues/", params; aws_config=aws_config)

"""
    list_alarms(alarm_model_name)
    list_alarms(alarm_model_name, params::Dict{String,<:Any})

Lists one or more alarms. The operation returns only the metadata associated with each
alarm.

# Arguments
- `alarm_model_name`: The name of the alarm model.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per request.
- `"nextToken"`: The token that you can use to return the next set of results.
"""
list_alarms(alarmModelName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/alarms/$(alarmModelName)"; aws_config=aws_config)
list_alarms(alarmModelName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/alarms/$(alarmModelName)", params; aws_config=aws_config)

"""
    list_detectors(detector_model_name)
    list_detectors(detector_model_name, params::Dict{String,<:Any})

Lists detectors (the instances of a detector model).

# Arguments
- `detector_model_name`: The name of the detector model whose detectors (instances) are
  listed.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per request.
- `"nextToken"`: The token that you can use to return the next set of results.
- `"stateName"`: A filter that limits results to those detectors (instances) in the given
  state.
"""
list_detectors(detectorModelName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)"; aws_config=aws_config)
list_detectors(detectorModelName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)", params; aws_config=aws_config)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetInstanceViewStatusesSummary <: SwaggerModel
    statusesSummary::Any # spec type: Union{ Nothing, Vector{VirtualMachineStatusCodeCount} } # spec name: statusesSummary

    function VirtualMachineScaleSetInstanceViewStatusesSummary(;statusesSummary=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetInstanceViewStatusesSummary, Symbol("statusesSummary"), statusesSummary)
        setfield!(o, Symbol("statusesSummary"), statusesSummary)
        o
    end
end # type VirtualMachineScaleSetInstanceViewStatusesSummary

const _property_map_VirtualMachineScaleSetInstanceViewStatusesSummary = Dict{Symbol,Symbol}(Symbol("statusesSummary")=>Symbol("statusesSummary"))
const _property_types_VirtualMachineScaleSetInstanceViewStatusesSummary = Dict{Symbol,String}(Symbol("statusesSummary")=>"Vector{VirtualMachineStatusCodeCount}")
Base.propertynames(::Type{ VirtualMachineScaleSetInstanceViewStatusesSummary }) = collect(keys(_property_map_VirtualMachineScaleSetInstanceViewStatusesSummary))
Swagger.property_type(::Type{ VirtualMachineScaleSetInstanceViewStatusesSummary }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetInstanceViewStatusesSummary[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetInstanceViewStatusesSummary }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetInstanceViewStatusesSummary[property_name]

function check_required(o::VirtualMachineScaleSetInstanceViewStatusesSummary)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetInstanceViewStatusesSummary }, name::Symbol, val)
end

ext: jl


text:
# syntax: proto3
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type Timestamp
    seconds::Int64
    nanos::Int32
    Timestamp(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type Timestamp
hash(v::Timestamp) = ProtoBuf.protohash(v)
isequal(v1::Timestamp, v2::Timestamp) = ProtoBuf.protoisequal(v1, v2)
==(v1::Timestamp, v2::Timestamp) = ProtoBuf.protoeq(v1, v2)

export Timestamp


text:
@testset "multiclass" begin
    y = [1, 2, 3, 1, 1, 2, 3]
    l = LogisticSoftMaxLikelihood(3)
    AGP.create_mapping!(l, y)
    @test sort(l.class_mapping) == [1, 2, 3]
    @test l.ind_mapping == Dict(1 => 1, 2 => 2, 3 => 3)
    Y = AGP.create_one_hot(l, y[1:3])
    @test Y == BitMatrix([1 0 0; 0 1 0; 0 0 1])

    l = LogisticSoftMaxLikelihood(2)
    @test_throws ErrorException AGP.create_mapping!(l, y)

    y = [1, 2, 1, 1]
    l = LogisticSoftMaxLikelihood(3)
    AGP.create_mapping!(l, y)
    @test l.class_mapping == [1, 2, 3]
    @test l.ind_mapping == Dict(1 => 1, 2 => 2, 3 => 3)
    Y = AGP.create_one_hot(l, y)
    @test Y == BitMatrix([1 0 0; 0 1 0; 1 0 0; 1 0 0])
    @test AGP.n_latent(l) == 3

    y = ["b", "a", "c", "a", "a"]
    l = LogisticSoftMaxLikelihood(3)
    AGP.create_mapping!(l, y)
    @test l.class_mapping == ["b", "a", "c"]
    @test l.ind_mapping == Dict("b" => 1, "a" => 2, "c" => 3)
    Y = AGP.create_one_hot(l, y)
    @test Y == BitMatrix([1 0 0; 0 1 0; 0 0 1; 0 1 0; 0 1 0])

    l = LogisticSoftMaxLikelihood(["a", "b", "c"])
    @test l.ind_mapping == Dict("a" => 1, "b" => 2, "c" => 3)
    Y = AGP.create_one_hot(l, y)
    @test Y == BitMatrix([0 1 0; 1 0 0; 0 0 1; 1 0 0; 1 0 0])

    l = LogisticSoftMaxLikelihood(3)
    Y = AGP.treat_labels!(y, l)
    @test l.class_mapping == ["b", "a", "c"]
    @test l.ind_mapping == Dict("b" => 1, "a" => 2, "c" => 3)
    @test Y == BitMatrix([1 0 0; 0 1 0; 0 0 1; 0 1 0; 0 1 0])
end


text:
"""
# SPNetworks.jl

A package for manipulating Sum-Product Networks.

## Introduction

A Sum-Product Network (SPN) is either

  - a univariate distribution (leaf node), or
  - a convex combination of SPNs (sum node), or
  - a product of SPNs.

SPNs are represented as rooted graphs; internal nodes represent convex
combinations and products; leaves represent distributions.

"""
module SPNetworks

# required packages

import Random

# Implements basic data types and helper functions
include("DataTypes/Nodes.jl")
include("DataTypes/SumProductNetworks.jl")

# Implements helper functions not dependent on SPNs
include("Utils.jl")

# For evaluating and sampling from SPN
include("Inference/Evaluation.jl")

# MAP Inference
include("Inference/MAP/MAP.jl")

# I/O functions
include("IO.jl")

# For learning the parameters of SPNs
include("Learning/ParameterLearning.jl")

# For learning the structure of SPNs
# include("Learning/StructureLearning.jl")
include("DataTypes/BayesianNets.jl")

# For generating Dense SPNs (for use e.g. in images)
# include("Learning/DenseSPNs.jl")

end # end of module


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Authors: Rafael de Pelegrini Soares
#*          Andrey Copat, Estefane S. Horn, Marcos L. Alencastro
#* $Id$
#*--------------------------------------------------------------------
#Needs to be reformulated
type HidraulicTurbine
	HidraulicTurbine()=begin
		NComp=outers.NComp
		PP=outers.PP
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components",
				:Lower=>1
			)),
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Molar Weight"
			)),(NComp)),
			efficiency(Dict{Symbol,Any}(
				:Brief=>"Turbine efficiency"
			)),
			efficiency(Dict{Symbol,Any}(
				:Brief=>"Brake efficiency"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Volumetric expansivity",
				:Unit=>"1/K"
			)),
			head(Dict{Symbol,Any}(
				:Brief=>"Head Developed"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Fluid Power"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Brake Power"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pressure Ratio"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			molweight(Dict{Symbol,Any}(
				:Brief=>"Mixture Molar Weight"
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Specific Mass"
			)),
			cp_mol(Dict{Symbol,Any}(
				:Brief=>"Heat Capacity"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0.05,
				:PosY=>0.0,
				:Symbol=>"_{in}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>0.65,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Work Outlet",
				:PosX=>1,
				:PosY=>0.46
			)),
			[
				:(Mwm = sum(Mw*Inlet.z)),
				:(rho = PP.LiquidDensity(Inlet.T,Inlet.P,Inlet.z)),
				:(Outlet.v = PP.VapourFraction(Outlet.T, Outlet.P, Outlet.z)),
				:(Cp = PP.LiquidCp(Inlet.T,Inlet.P,Inlet.z)),
				:(Outlet.P = Inlet.P * Pratio),
				:(Outlet.P = Inlet.P - Pdrop),
				:(FPower * rho = -Pdrop * Inlet.F * Mwm),
				:(BPower = FPower * Eff),
				:(BPower = WorkOut),
				:((Outlet.T - Inlet.T) * rho * Cp = (Outlet.h - Inlet.h) * rho + Pdrop * Mwm * (1-Beta*Inlet.T)),
				:((Outlet.h - Inlet.h) * rho = -Pdrop * Mwm),
				:(Outlet.F = Inlet.F),
				:(Outlet.z = Inlet.z),
				:(Head = Outlet.h - Inlet.h),
			],
			[
				"Calculate Mwm for Inlet Mixture","Calculate rho using a External Physical Properties Routine","Calculate Outlet Vapour Fraction","Calculate Cp Using a External Physical Properties Routine","Pressure Ratio","Pressure Drop","Calculate Fluid Power","Calculate Brake Power","","Calculate Outlet Temperature","Calculate Outlet Enthalpy","Molar Balance","Calculate Outlet Composition","Calculate Head",
			],
			[:NComp,:PP,:Mw,],
			[:Eff,:Meff,:Beta,:Head,:FPower,:BPower,:Pratio,:Pdrop,:Mwm,:rho,:Cp,:Inlet,:Outlet,:WorkOut,]
		)
	end
	NComp::DanaInteger
	PP::DanaPlugin
	Mw::Array{molweight}
	Eff::efficiency
	Meff::efficiency
	Beta::positive
	Head::head
	FPower::power
	BPower::power
	Pratio::positive
	Pdrop::press_delta
	Mwm::molweight
	rho::dens_mass
	Cp::cp_mol
	Inlet::stream
	Outlet::stream
	WorkOut::power
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export HidraulicTurbine
function set(in::HidraulicTurbine)
	Mw = PP.MolecularWeight()

end
function setEquationFlow(in::HidraulicTurbine)
	#Mixtures Properties
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
end
function atributes(in::HidraulicTurbine,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/HidraulicTurbine"
	fields[:Brief]="Testing Model of a Hidraulic Turbine."
	drive!(fields,_)
	return fields
end
HidraulicTurbine(_::Dict{Symbol,Any})=begin
	newModel=HidraulicTurbine()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# Inside make.jl
push!(LOAD_PATH,"../src/")
using BlockDates
using Documenter

makedocs(
         sitename = "BlockDates.jl",
         modules  = [BlockDates],
         pages=[
                "Home" => "index.md"
               ])
deploydocs(;
    repo="github.com/EconometricsBySimulation/BlockDates.jl",
)


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("MLJTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# # Simple example of a homogeneous ensemble using learning networks
# In this simple example, no bagging is used, so every atomic model# gets the same learned parameters, unless the atomic model training# algorithm has randomness, eg, DecisionTree with random subsampling# of features at nodes.
# Note that MLJ has a built in model wrapper called `EnsembleModel`# for creating bagged ensembles with a few lines of code.
# ## Definition of composite model type
using MLJ, PyPlot
import Statistics

# learning network (composite model spec):
Xs = source()
ys = source(kind=:target)

atom = @load DecisionTreeRegressor
atom.n_subfeatures = 4 # to ensure diversity among trained atomic models

machines = (machine(atom, Xs, ys) for i in 1:100)

# overload `mean` for nodes:
Statistics.mean(v...) = mean(v)
Statistics.mean(v::AbstractVector{<:AbstractNode}) = node(mean, v...)

yhat = mean([predict(m, Xs) for  m in machines]);

# new composite model type and instance:
one_hundred_models = @from_network OneHundredModels(atom=atom) <= yhat

# ## Application to data
X, y = @load_boston;

# tune regularization parameter for a *single* tree:
r = range(atom,
          :min_samples_split,
          lower=2,
          upper=100, scale=:log)

mach = machine(atom, X, y)

curve = learning_curve!(mach,
                        range=r,
                        measure=mav,
                        resampling=CV(nfolds=9),
                        verbosity=0)

plot(curve.parameter_values, curve.measurements)
xlabel(curve.parameter_name)

# \fig{e1.svg}
# tune regularization parameter for all trees in ensemble simultaneously:
r = range(one_hundred_models,
          :(atom.min_samples_split),
          lower=2,
          upper=100, scale=:log)

mach = machine(one_hundred_models, X, y)

curve = learning_curve!(mach,
                        range=r,
                        measure=mav,
                        resampling=CV(nfolds=9),
                        verbosity=0)

plot(curve.parameter_values, curve.measurements)
xlabel(curve.parameter_name)

# \fig{e2}

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
# License for this file: MIT (expat)
# Copyright 2017-2019, DLR Institute of System Dynamics and Control

module Modia3D

const path = dirname(dirname(@__FILE__))   # Absolute path of package directory
const Version = "0.4.1-dev"
const Date = "2019-11-12"

println("\nImporting Modia3D Version $Version ($Date)")

# Abstract types
#abstract type AbstractComponentWithVariables end                    # Component that has potentially Variables

abstract type AbstractFlange end              # Flanges are associated with joints, frames, ... (are for connectin parts)

abstract type AbstractObject3Ddata end                                 # Data associated with one Object3D
abstract type AbstractVisualElement  <: AbstractObject3Ddata     end   # Visual element associated with one Object3D
abstract type AbstractGeometry       <: AbstractVisualElement end   # Geometry type
abstract type AbstractCamera         <: AbstractVisualElement end   # Camera type
abstract type AbstractHeadUpDisplay  <: AbstractVisualElement end   # Head-Up-Display type
abstract type AbstractPathDefinition <: AbstractVisualElement end   # Path definition type
abstract type AbstractEffect         <: AbstractVisualElement end   # Effect type

abstract type AbstractSolidGeometry     end   # Immutable shape type that  has a volume and optionally mass. Can be used in collisions.
abstract type AbstractSolidMaterial     end   # Material properties of a solid (e.g. density)
abstract type AbstractContactMaterial   end   # Contact properties of a solid (e.g. spring constant)
abstract type AbstractContactDetection  end   # Contact detection type
abstract type AbstractGravityField      end   # Gravity field type
abstract type AbstractRenderer          end   # Renderer type
abstract type AbstractDLR_VisualizationRenderer <: AbstractRenderer end   # Community or Professional edition of DLR_Visualization renderer
#abstract type AbstractSignal            end   # Signals
#abstract type AbstractBus               end   # Bus: is a collection of signals

abstract type AbstractMassProperties    end

abstract type AbstractKeys    end
abstract type AbstractContactPairMaterial end # Constants needed to compute the contact response between two objects

import ModiaMath
abstract type AbstractAssemblyInternal  <: ModiaMath.AbstractComponentInternal end
abstract type AbstractAssemblyComponent <: ModiaMath.AbstractComponentWithVariables end # Objects that can be stored in an Assembly (has elements name and within)
abstract type AbstractAssembly          <: AbstractAssemblyComponent end  # Assembly object (has elements name, within, referenceObject3D and scene)
abstract type AbstractTwoObject3DObject <: AbstractAssemblyComponent end  # Object related to two Object3Ds
abstract type AbstractJoint             <: AbstractTwoObject3DObject end  # Constraint between two Object3Ds
abstract type AbstractRevolute          <: AbstractJoint             end  # Revolute joint
abstract type AbstractPrismatic         <: AbstractJoint             end  # Prismatic joint

abstract type AbstractForceTorque <: AbstractAssemblyComponent end  # Object computing the torque of a revolute joint

abstract type AbstractSignal <: AbstractAssemblyComponent end # It is a Signal with one output variable
abstract type AbstractSignalAdaptor  <: AbstractAssemblyComponent end
abstract type AbstractForceAdaptor   <: AbstractAssemblyComponent end

# Enumerations
@enum Ternary      True False Inherited
@enum AnalysisType         KinematicAnalysis QuasiStaticAnalysis           DynamicAnalysis
@enum VariableAnalysisType AllAnalysis       QuasiStaticAndDynamicAnalysis OnlyDynamicAnalysis NotUsedInAnalysis

# Used renderer (actual value is defined with __init__() below)
const renderer = Vector{AbstractRenderer}(undef,1)

# Include sub-modules
include(joinpath("Basics"          , "_module.jl"))
include(joinpath("Signals"         , "_module.jl"))
include(joinpath("ForceElements"   , "_module.jl"))
include(joinpath("Graphics"        , "_module.jl"))
include(joinpath("Solids"          , "_module.jl"))
include(joinpath("Composition"     , "_module.jl"))
include(joinpath("renderer"        , "DLR_Visualization"  , "_module.jl"))
include(joinpath("renderer"        , "NoRenderer"         , "_module.jl"))
include(joinpath("contactDetection", "ContactDetectionMPR", "_module.jl"))

# Make symbols available that have been exported in sub-modules
using  .Basics
using  .Signals
using  .ForceElements
using  .Solids
using  .Graphics
using  .Composition
import .DLR_Visualization
import .NoRenderer
const Pipe     = Graphics.Pipe        # Pipe cannot be directly exported, due to a conflict with Base.Pipe
const connect  = Composition.connect  # connect cannot be directly exported, due to a conflict with Base.connect

# Called implicitely at the first import/using of Modia3D (when loading Modia3D to the current Julia session)
function __init__()
    info = DLR_Visualization.getSimVisInfo()
    (directory, dll_name, isProfessionalEdition, isNoRenderer) = info

    if isNoRenderer
        renderer[1] = NoRenderer.DummyRenderer(info)
    elseif isProfessionalEdition
        renderer[1] = DLR_Visualization.ProfessionalEdition(info)
    else
        renderer[1] = DLR_Visualization.CommunityEdition(info)
    end
end

export @assembly
export @signal
export @bus
export @forceElement
export Object3D

export Sphere, Ellipsoid, Box, Cylinder, Capsule, Beam, Cone, Pipe
export Spring, GearWheel, CoordinateSystem, Grid, FileMesh

export Solid, SolidFileMesh, SolidWithConvexDecomposition
export SolidSphere, SolidEllipsoid, SolidBox, SolidCylinder, SolidCapsule, SolidBeam, SolidCone, SolidPipe
export MassProperties
export Fixed, Revolute, Prismatic

export UniformGravityField, PointGravityField, NoGravityField
export Material
export SceneOptions, SimulationModel
export visualizeAssembly!, print_ModelVariables

export PTP_path, getPosition!, getPosition, getIndex, plotPath

# Add import clauses used in examples and test
import StaticArrays
import Unitful
import ModiaMath
import LinearAlgebra
import Test

end # module


text:
"""
    Sherlock(optimizer, ϵ::Float64)

Sherlock combines local and global search to estimate the range of the output node.

# Problem requirement
1. Network: any depth, ReLU activation, single output
2. Input: hpolytope and hyperrectangle
3. Output: hyperrectangle (1d interval)

# Return
`CounterExampleResult` or `ReachabilityResult`

# Method
Local search: solve a linear program to find local optima on a line segment of the piece-wise linear network.
Global search: solve a feasibilty problem using MILP encoding (default calling NSVerify).
- `optimizer` default `GLPKSolverMIP()`
- `ϵ` is the margin for global search, default `0.1`.

# Property
Sound but not complete.

# Reference
[S. Dutta, S. Jha, S. Sanakaranarayanan, and A. Tiwari,
"Output Range Analysis for Deep Neural Networks,"
*ArXiv Preprint ArXiv:1709.09130*, 2017.](https://arxiv.org/abs/1709.09130)

[https://github.com/souradeep-111/sherlock](https://github.com/souradeep-111/sherlock)
"""
@with_kw struct Sherlock
    optimizer::AbstractMathProgSolver = GLPKSolverMIP()
    ϵ::Float64                        = 0.1
end

function solve(solver::Sherlock, problem::Problem)
    (x_u, u) = output_bound(solver, problem, :max)
    (x_l, l) = output_bound(solver, problem, :min)
    bound = Hyperrectangle(low = [l], high = [u])
    reach = Hyperrectangle(low = [l - solver.ϵ], high = [u + solver.ϵ])
    return interpret_result(reach, bound, problem.output, x_l, x_u) # This function is defined in bab.jl
end

function output_bound(solver::Sherlock, problem::Problem, type::Symbol)
    opt = solver.optimizer
    x = sample(problem.input)
    while true
        (x, bound) = local_search(problem, x, opt, type)
        bound += ifelse(type == :max, solver.ϵ, -solver.ϵ)
        (x_new, bound_new, feasibile) = global_search(problem, bound, opt, type)
        feasibile || return (x, bound)
        (x, bound) = (x_new, bound_new)
    end
end

# Choose the first vertex
function sample(set::AbstractPolytope)
    x = vertices_list(set)
    return x[1]
end

function local_search(problem::Problem, x::Vector{Float64}, optimizer::AbstractMathProgSolver, type::Symbol)
    nnet = problem.network
    act_pattern = get_activation(nnet, x)
    gradient = get_gradient(nnet, x)
    model = Model(solver = optimizer)
    neurons = init_neurons(model, nnet)
    add_set_constraint!(model, problem.input, first(neurons))
    encode_lp!(model, nnet, neurons, act_pattern)
    o = gradient * neurons[1]
    index = ifelse(type == :max, 1, -1)
    @objective(model, Max, index * o[1])
    solve(model, suppress_warnings = true)
    x_new = getvalue(neurons[1])
    bound_new = compute_output(nnet, x_new)
    return (x_new, bound_new[1])
end

function global_search(problem::Problem, bound::Float64, optimizer::AbstractMathProgSolver, type::Symbol)
    index = ifelse(type == :max, 1.0, -1.0)
    h = HalfSpace([index], index * bound)
    output_set = HPolytope([h])
    problem_new = Problem(problem.network, problem.input, output_set)
    solver  = NSVerify(optimizer = optimizer)
    result  = solve(solver, problem_new)
    if result.status == :UNSAT
        x = result.counter_example
        bound = compute_output(problem.network, x)
        return (x, bound[1], true)
    else
        return ([], 0.0, false)
    end
end


text:
using Morbit
using Test
using Logging

@testset "Save & Load AbstractConfig" begin

    # AbstractConfig
    empty_cfg = Morbit.DefaultConfig();
    fn1 = string( tempname(), ".jld2" )
    pt1 = save_config( fn1, empty_cfg )

    @test fn1 == pt1
    @test isfile(pt1)

    loaded_empty_cfg = load_config( fn1 )
    @test empty_cfg == loaded_empty_cfg

    algo_cfg = AlgorithmConfig();
    fn2 = tempname()
    pt2 = save_config( fn2, algo_cfg )

    @test string(fn2, ".jld2") == pt2
    @test isfile(pt2)

    loaded_algo_config = load_config( pt2 )
    CFG_SAME = true
    for field in fieldnames( AlgorithmConfig )
        if getfield( algo_cfg, field ) ≠ getfield( loaded_algo_config, field )
            CFG_SAME = false
        end
    end
    @test CFG_SAME

    # try to load non-existent file
    # use NullLogger to suppress threatening messages
    with_logger(NullLogger()) do
        @test isnothing(load_config(fn2))
    end
end

@testset "Save & Load Database" begin

    # AbstractConfig
    empty_db = Morbit.NoDB();
    fn1 = string( tempname(), ".jld2" )
    pt1 = save_database( fn1, empty_db )

    @test fn1 == pt1
    @test isfile(pt1)

    loaded_empty_db = load_database( pt1 )
    @test empty_db == loaded_empty_db

    fn2 = tempname()
    nonempty_db = Morbit.ArrayDB();
    Morbit.add_result!(nonempty_db,
        Morbit.init_res(
            Morbit.Result,
            rand(2),
            rand(3)
        )
    )
    pt2 = save_database(fn2, nonempty_db)
    @test string(fn2, ".jld2") == pt2
    @test isfile(pt2)

    loaded_nonempty_db = load_database( pt2 )

    DB_SAME = true
    for field in fieldnames( Morbit.ArrayDB )
        if getfield( nonempty_db, field ) ≠ getfield( loaded_nonempty_db, field )
            @show getfield(nonempty_db, field)
            @show getfield(loaded_nonempty_db, field)
            DB_SAME = false
            break;
        end
    end
    @test DB_SAME
end

# TODO test save / load of AbstractIterData


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

using FiniteStateTransducers
using Test, Random
using DataStructures
Random.seed!(123)

@testset "FiniteStateTransducers.jl" begin
  @testset "Semirings" begin
    include("semirings.jl")
  end
  @testset "FiniteStateTransducers constructors and utils" begin
    include("wfst.jl")
    include("properties.jl")
    include("io_wfst.jl")
  end
  @testset "algorithms" begin
    include("paths.jl") #Paths, BFS, paths
    include("dfs.jl") #rm_state, DFS, get_scc, connect, rm_eps
    include("shortest_distance.jl")
    include("topsort.jl")
    include("rm_eps.jl")
    include("closure.jl")
    include("inv.jl")
    include("proj.jl")
    include("reverse.jl")
    include("compose.jl")
    include("cat.jl")
    include("determinize.jl")
    include("wfst2tr.jl")
  end
end

#@testset "FiniteStateTransducers.jl" begin
#end


text:
using BSplineKit
using BSplineKit.BSplines: multiplicity

using BandedMatrices
using LinearAlgebra
using Test

eval_poly(x::Number, P) = @evalpoly(x, P...)

# Test a polynomial of degree k - 1.
# The splines should approximate the polynomial (and its derivatives) perfectly.
# This is also used to test the SplineInterpolations module.
function test_polynomial(x, ::BSplineOrder{k}) where {k}
    # Coefficients of polynomial of degree k - 1 (see ?evalpoly).
    # P(x) = -1 + 2x - 3x^2 + 4x^3 - ...
    P = ntuple(d -> (-d)^d, Val(k))
    P′ = ntuple(d -> d * P[d + 1], Val(k - 1))    # derivative
    Pint = (0, ntuple(d -> P[d] / d, Val(k))...)  # antiderivative

    # Interpolate polynomial at `x` locations.
    y = eval_poly.(x, Ref(P))
    itp = @inferred interpolate(x, y, BSplineOrder(k))

    @test startswith(repr(itp), "SplineInterpolation containing")

    S = spline(itp)
    @test length(S) == length(x)

    @testset "Interpolations" begin
        @test itp.(x) ≈ y

        let x = (x[2] + x[3]) / 3
            @test itp(x) == S(x)  # these are equivalent
        end

        @test order(itp) === order(S)
        @test basis(itp) === basis(S)
        @test knots(itp) === knots(S)
        @test coefficients(itp) === coefficients(S)
        @test integral(itp) == integral(S)
        @test diff(itp, Derivative(1)) == diff(S, Derivative(1)) == Derivative() * itp

        # "incompatible lengths of B-spline basis and collocation points"
        @test_throws(
            DimensionMismatch,
            SplineInterpolation(undef, basis(S), x[1:4], eltype(S)),
        )

        let B = basis(S)
            xs = collocation_points!(Vector{Float32}(undef, length(B)), B)
            itp = @inferred SplineInterpolation(undef, B, xs)
            @test eltype(itp) === Float32
        end

        # "input data has incorrect length"
        @test_throws DimensionMismatch interpolate!(itp, rand(length(x) - 1))
    end

    S′ = diff(S, Derivative(1))
    @test Derivative() * S == S′  # alternative notation
    Sint = integral(S)

    a, b = boundaries(basis(S))

    @test Sint(a) == 0  # this is an arbitrary choice
    Pint_a = @evalpoly(a, Pint...)

    # Compare values on a finer grid.
    let Nx = 9 * length(S) + 42
        x = LinRange(a, b, Nx)
        @test all(@evalpoly(x, P...) ≈ S(x) for x in x)
        @test all(@evalpoly(x, P′...) ≈ S′(x) for x in x)
        @test all(@evalpoly(x, Pint...) - Pint_a ≈ Sint(x) for x in x)
    end

    nothing
end

function test_splines(B::BSplineBasis, knots_in)
    k = order(B)
    t = knots(B)

    @testset "Knots (k = $k)" begin
        let (ka, kb) = multiplicity.(Ref(t), (1, length(t)))
            @test ka == kb == k
        end

        @test @views all(t[1:k] .== knots_in[1])
        @test @views all(t[(end - k + 1):end] .== knots_in[end])
        @test @views t[(k + 1):(end - k)] == knots_in[2:(end - 1)]
    end

    @testset "B-splines (k = $k)" begin
        N = length(B)
        @test_throws DomainError evaluate(B, 0, 0.2)
        @test_throws DomainError evaluate(B, N + 1, 0.2)

        @test_throws BoundsError B[0]
        @test_throws BoundsError B[N + 1]

        # In segment [t[k + 1], t[k + 2]], non-zero B-splines are 2:(k + 1).
        @test @inferred(nonzero_in_segment(B, k + 1)) == 2:(k + 1)

        # Fewer B-splines near the boundaries. This is usually not a problem
        # because knots are repeated there...
        @test nonzero_in_segment(B, 1) == 1:1
        @test nonzero_in_segment(B, N) == (N - k + 1):N
        @test nonzero_in_segment(B, N + 1) == (N - k + 2):N
        @test nonzero_in_segment(B, N + k - 1) == N:N
        @test isempty(nonzero_in_segment(B, N + k))

        # These intervals are not defined
        @test isempty(nonzero_in_segment(B, 0))
        @test isempty(nonzero_in_segment(B, length(knots(B))))

        # Verify values at the boundaries.
        a, b = boundaries(B)
        @test evaluate(B, 1, a) == 1.0
        @test evaluate(B, N, b) == 1.0

        @test startswith(
            repr(B),
            "$N-element BSplineBasis of order $k",
        )

        xs = range(a, b; length = 4N + 1)
        us = similar(xs)
        vs = similar(us)
        i = N >> 2
        evaluate!(us, B, i, xs)
        @test us == evaluate(B, i, xs)
        @test us == BasisFunction(B, i).(xs)
        @test us == B[i].(xs)
        @test us == map!(B[i], vs, xs)
        @test us ≈ B[i, Float32].(xs)

        @testset "iterate" begin
            for (i, b) in enumerate(B)
                @test b === B[i]
            end
        end
    end

    xcol = collocation_points(B, method=Collocation.AvgKnots())
    @test xcol[1] == knots_in[1]
    @test xcol[end] == knots_in[end]

    C = @inferred collocation_matrix(B, xcol)

    let N = size(C, 1)
        T = eltype(C)
        l, u = bandwidths(C)
        @test startswith(
            sprint(show, MIME("text/plain"), C),
            "$N×$N CollocationMatrix{$T} with bandwidths ($l, $u):\n",
        )
    end

    @testset "LU factorisation" begin
        T = eltype(C)
        @test T === Float64
        @test C isa CollocationMatrix{T}
        F = @inferred factorize(C)
        @test F isa LU{T, <:CollocationMatrix{T}}
        @test F.L * F.U ≈ C
        @test F.P == I  # permutation matrix = identity matrix (no pivoting)
        @test F.p == 1:size(C, 1)  # no row permutations

        y = sin.(xcol)
        u = F \ y
        @test C * u ≈ y

        let v = similar(u, length(u) - 2)
            @test_throws DimensionMismatch ldiv!(v, F, y)
        end
    end

    @testset "Spline (k = $k)" begin
        @testset "Polynomial" begin
            test_polynomial(xcol, BSplineOrder(k))
        end

        # Generate data at collocation points and get B-spline coefficients.
        ucol = cos.(xcol)
        coefs = C \ ucol
        @test C * coefs ≈ ucol

        @inferred Spline(B, coefs)
        S = Spline(B, coefs)
        @test length(S) == length(B)
        @test all(S.(xcol) .≈ ucol)
        @test coefficients(S) === coefs
        @test diff(S, Derivative(0)) === S
        @test_nowarn show(devnull, S)
        @test Splines.parent_spline(S) === S

        # Broadcasting
        f(S, x) = S(x)
        @test f.(S, xcol) == S.(xcol)

        # Create new spline, then compare it to S.
        let P = Spline(undef, B)
            cp = coefficients(P)
            fill!(cp, 0)
            @test P != S
            @test !(P ≈ S)

            copy!(cp, coefs)  # copy coefficients of S
            @test P == S
            @test P ≈ S

            @test copy(S) == P
        end
    end

    nothing
end

function test_splines(::BSplineOrder{k}) where {k}
    gauss_lobatto_points(N) = [-cos(π * n / N) for n = 0:N]
    breaks = (
        -1:0.05:1,
        gauss_lobatto_points(10 + k),
    )
    for x in breaks
        @inferred BSplineBasis(BSplineOrder(k), copy(x))
        @inferred BSplineBasis(BSplineOrder(k), copy(x), augment=Val(false))
        @inferred BSplineBasis(BSplineOrder(k), copy(x), augment=Val(true))
        @inferred (() -> BSplineBasis(k, copy(x)))()
        g = BSplineBasis(k, copy(x))

        @testset "BSplineBasis equality" begin
            h = BSplineBasis(k, copy(x))
            @test g == h

            h = BSplineBasis(BSplineOrder(k + 1), copy(x))
            @test g != h

            h = BSplineBasis(k, x .+ 1)
            @test g != h
        end

        @test order(g) == k
        test_splines(g, x)
    end
    nothing
end

@testset "Splines" begin
    test_splines(BSplineOrder(4))
    test_splines(BSplineOrder(5))
end


text:
function solveBooleanReg(data_dictionary, initial_state, initial_input1, initial_input2)

#bring in the data that I need to evaluate the regulatory states
regulatoryRules = data_dictionary["regulatoryRules"]
regulatoryInputs1 = data_dictionary["regulatoryInputs1"]
regulatoryInputs2 = data_dictionary["regulatoryInputs2"]
flux_bounds = data_dictionary["flux_bounds_array"]
regulatoryGenes = data_dictionary["regulatoryGenes"]
genes = data_dictionary["list_of_gene_symbols"]
flux_names = data_dictionary["list_of_reaction_name_strings"]
final_input1 = false.*ones(length(data_dictionary["regulatoryInputs1"])); #initialize
s = data_dictionary["stoichiometric_matrix"]

initial_state = initial_state .== 1
initial_input1 = initial_input1 .== 1
initial_input2 = initial_input2 .== 1
initial_regstate = vcat(initial_state, initial_input1, initial_input2) #Full list of states

#This determines external metabolite levels using flux bounds (input 1)
for i = 1:length(regulatoryInputs1)
    flux_n = regulatoryInputs1[i]; #get the metabolite in the first regulatory input
    metID = findfirst(isequal(flux_n), data_dictionary["list_of_metabolite_symbols"]) #find the index of that metabolite
    potfluxID = findall(isequal(~0), s[metID[1],:]) #find the fluxes associated with that metabolite
#so now we only need to select actual exchanges, meaning only one -1 or 1.
    for j = 1:length(potfluxID)
        if abs(sum(s[:,potfluxID[j]])) == 1
            rxnID = potfluxID[j]
            #println(rxnID)
                if flux_bounds[rxnID, 1] < 0 #if exchange is allowed
                    final_input1[i] = true;
                    #println(flux_bounds[rxnID,1])
                    #println("True")
                else
                    final_input1[i] = false; #if exchange is not allowed
                    #println(flux_bounds[rxnID,1])
                    #println("false")
                end
        end
    end

end #now we have new input1 states

#apply initial regulatory state to model and obtain FBA
drGenes = []
for i = 1:length(regulatoryGenes)
    if initial_state[i] == false
        push!(drGenes, regulatoryGenes[i])
    end
end
drGenes = intersect(genes,drGenes)
data_dictionary_deleted = deleteModelGenes(data_dictionary, drGenes)[1]
#data_dictionary_deleted = optimize_specific_growth_rate(data_dictionary_deleted) #placeholder for now
(objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag, status_flag) = calculate_optimal_flux_distribution(data_dictionary_deleted); #change this to 2 version if standalone
#println("Succ:", calculated_flux_array[87])
#Now, we need to check regulatory inputs that are dependent on fluxes
final_input2 = true.*ones(length(data_dictionary["regulatoryInputs2"]))
#println(calculated_flux_array)
fluxInd = indexin(regulatoryInputs2, flux_names) #figure out where these inputs are in the list of fluxes
for i = 1:length(regulatoryInputs2)
    if calculated_flux_array[fluxInd[i]] == 0 #if there was no flux in the regulatory input 2 flux
        final_input2[i] = false
    end
end

#Now that we have regulatory inputs 1 and 2 set, it's time to set the state of genes for next time

geneStates = calculate_full_rules_vector(data_dictionary, initial_regstate)
#geneStates = geneStates .== 1
final_state = geneStates #the final state of the regulatory system

return(final_state,final_input1,final_input2)
#so the next time around, the final_state will differ from the initial state as
#we now are inputting new inputs1 and inputs2.
end


text:
# small test
using LightGraphs; using D3GraphPlot; d3plot(Graph(4,6))


text:
using VLOptionsModeling
using Test

# -- Model creation tests ------------------------------------------------------ #
function run_default_test()
    return true
end
# ------------------------------------------------------------------------------- #

@testset "default_test_set" begin
    @test run_default_test() == true
end


text:
@testset "Verify Jacobian of off-diagonal entries" begin
    Nx = 4
    Ne = 100

    X = randn(Nx, Ne) .* randn(Nx, Ne)
    λ = 0.0
    δ = 0.0
    γ = 2.0

    for p in [0; 1; 2; 3]
        if p == 0
            coeff_off = randn((Nx-1)*(p+1))
            coeff_diag = rand(p+2)
        else
            coeff_off = randn((Nx-1)*(p+1))
            coeff_diag = rand(p+3)
        end

        coeff = vcat(coeff_off, coeff_diag)
        dJ = zero(coeff_off)
        order = p*ones(Int64, Nx)
        Cfull = SparseRadialMapComponent(Nx, order)
        center_std!(Cfull, sort(X; dims = 2); γ = γ)
        modifycoeff!(Cfull, coeff)
        ψ_off, ψ_diag, ψ_∂k = compute_weights(Cfull, X)

        J(x) = (1/(2*Ne))*norm((ψ_off*x) + ψ_diag*coeff_diag[2:end] .+ coeff_diag[1])^2

        dJautodiff = ForwardDiff.gradient(x-> J(x), coeff_off)

        rhs = zeros(Ne)
        cache = zeros(Ne)
        mul!(rhs, ψ_diag, coeff_diag[2:end])
        rhs .+= coeff_diag[1]
        rmul!(rhs, -1.0)
        gradient_off!(dJ, cache, ψ_off, coeff_off, rhs, Ne)

        @test isapprox(dJ, dJautodiff, atol = 1e-10)
    end
end

@testset "Test greedyfit for a fixed number of family of features" begin

    for Nx=1:5
        Ne = 100

        X = randn(Nx, Ne) .* randn(Nx, Ne)

        λ = 0.0
        δ = 1e-8
        γ = 2.0

        for pdiag=0:3
            for poff=0:3
                for k=0:Nx-1
                    Cgreedy, ϵgreedy = greedyfit(Nx, poff, pdiag, deepcopy(X), k, λ, δ, γ)
                    C = SparseRadialMapComponent(Nx, Cgreedy.p)
                    center_std!(C, sort(deepcopy(X); dims = 2); γ = γ)
                    C, error = optimize(C, X, nothing, λ, δ)
                    @test isapprox(ϵgreedy[end], negative_likelihood(C, X, λ, δ), atol = 1e-10)
                    @test norm(Cgreedy.coeff - C.coeff)<1e-10
                end
            end
        end

        for poff=0:3
            for pdiag=0:3
                for k=0:Nx-1
                    Cgreedy, ϵgreedy = greedyfit(Nx, poff, pdiag, deepcopy(X), k, λ, δ, γ)
                    C = SparseRadialMapComponent(Nx, Cgreedy.p)
                    center_std!(C, sort(deepcopy(X); dims = 2); γ = γ)
                    C, error = optimize(C, X, nothing, λ, δ)
                    @test isapprox(ϵgreedy[end], negative_likelihood(C, X, λ, δ), atol = 1e-10)
                    @test norm(Cgreedy.coeff - C.coeff)<1e-10
                end
            end
        end
    end
end


text:
using SpecialFunctions
const base_unary_real = (
    (:Base, :acos, (-1, 1)),
    (:Base, :asin, (-1, 1)),
    (:Base, :atan, (-50, 50)),
    (:Base, :cos, (-1000, 1000)),
    (:Base, :sin, (-1000, 1000)),
    (:Base, :tan, (-1000, 1000)),
    (:Base, :acosh, (1, 1000)),
    (:Base, :asinh, (-1000, 1000)),
    (:Base, :atanh, (-1, 1)),
    (:Base, :cosh, (0, 89.415985f0)),
    (:Base, :sinh, (-89.415985f0, 89.415985f0)),
    (:Base, :tanh, (-8.66434f0, 8.66434f0)),
    (:Base, :cbrt, (-1000, 1000)),
    (:Base, :sqrt, (0, 1000)),
    (:Base, :exp, (-88.72284f0, 88.72284f0)),
    (:Base, :expm1, (-88.72284f0, 88.72284f0)),
    (:Base, :log, (0, 1000)),
    (:Base, :log10, (0, 1000)),
    (:Base, :log1p, (-1, 1000)),
    (:Base, :abs, (-1000, 1000)),
    (:Base, :abs2, (-1000, 1000)),
    (:Base, :ceil, (-1000, 1000)),
    (:Base, :floor, (-1000, 1000)),
    (:Base, :round, (-1000, 1000)),
    (:Base, :trunc, (-1000, 1000)),
    (:Base, :cis, (-1000, 1000)),
    (:SpecialFunctions, :erf, (-3.8325067f0, 3.8325067f0)),
    (:SpecialFunctions, :erfc, (-3.7439213f0, 10.019834f0)),
    (:SpecialFunctions, :erfinv, (-1, 1)),
    (:SpecialFunctions, :erfcinv, (0, 2)),
    (:SpecialFunctions, :lgamma, (0, 1000)),
    (:SpecialFunctions, :gamma, (0, 36))
)

const base_binary_real = (
    (:Base, :atan, (-1, 1), (-1, 1)),
    (:Base, :hypot, (-1000, 1000), (-1000, 1000)),
    # (getfield(Base, :./), (-1000, 1000), (-1000, 1000)),
    # (getfield(Base, :.^), (0, 100), (-5, 20))
)

const base_unary_complex = (
    (:Base, :acos, (-1, 1)),
    (:Base, :asin, (-1, 1)),
    (:Base, :acosh, (1, 1000)),
    (:Base, :asinh, (-1000, 1000)),
    (:Base, :sqrt, (0, 1000)),
    (:Base, :exp, (-88.72284f0, 88.72284f0)),
    (:Base, :log, (0, 1000)),
    (:Base, :abs, (-1000, 1000)),
    (:Base, :angle, (-1000, 1000)),
    (:Base, :conj, (-1000, 1000)),
    # (atan, (-50, 50)),
    # (cos, (-10, 10)),
    # (sin, (-10, 10)),
    # (tan, (-10, 10)),
    # (atanh, (-1, 1)),
    # (cosh, (0, 89.415985f0)),
    # (sinh, (-89.415985f0, 89.415985f0)),
    # (tanh, (-8.66434f0, 8.66434f0)),
    # (log10, (0, 1000)),
    # (cis, (-1000, 1000))
)

# const base_binary_complex = (
#     # (getfield(Base, :./), (-1000, 1000), (-1000, 1000)),
#     # ((:.^), (0, 100), (-2, 10))
# )

function randindomain(t::Type{T}, n, domain) where {T<:Real}
    d1 = convert(t, domain[1])
    d2 = convert(t, domain[2])
    ddiff = d2 - d1
    @assert isfinite(ddiff)
    v = rand(t, n)
    for i = 1:length(v)
        v[i] = v[i]*ddiff+d1
    end
    v
end

function randindomain(t::Type{T}, n, domain) where {T<:Complex}
    d1 = convert(t, domain[1])
    d2 = convert(t, domain[2])
    ddiff = d2 - d1
    @assert isfinite(ddiff)
    v = rand(t, 2*n)
    for i = 1:length(v)
        v[i] = v[i]*ddiff+d1
    end
    v
    # reinterpret(t, v)
end


text:
module DataStructuresAlgorithms

using Random

export Stack, DynamicArray
export linearsearch, binarysearch
export selectionsort!, insertionsort!, mergesort, quicksort!
export SLinkedList, DLinkedList, addlast!, addfirst!, add!, removelast!, removefirst!, remove!, peeklast, peekfirst
export Queue, enqueue!, dequeue!, peek, peekfirst, peeklast
export RMHeap, BNode, add!, remove!, meld!, MinPQ, add!, isheap, greater, getmin
export BST, BSTNode, printsort, isbst, maximum_node
export SCHash, SCLHash, OAHash, compress

import Base.push!, Base.pop!, Base.isempty, Base.peek, Base.show
import Base.getindex, Base.resize!, Base.setindex!, Base.:(==)
import Base.iterate, Base.eltype, Base.length, Base.size, Base.get!, Base.put!
import Base.delete!, Base.collect, Base.minimum, Base.maximum, Base.get, Base.hash

include("dynamicarray.jl")
include("stack.jl")
include("list.jl")
include("queue.jl")
include("sort.jl")
include("search.jl")
include("search.jl")
include("heap.jl")
include("tree.jl")
include("hash.jl")

end # module


text:
"""
# Coordinate conversion with the Gauss-Krüger Projection
CoordinateConverterGK.jl is a Julia package for conversion between
the Cartesian coordinates and the geographic coordinates for a point with
the Gauss-Krüger Projection.
"""
module CoordinateConverterGK
    include("defconstant.jl")
    include("arclength.jl")
    include("yx2latlon.jl")
    include("latlon2yx.jl")
    export yx2latlon, xy2lonlat
    export yx2latlon_ja, xy2lonlat_ja
    export latlon2yx, lonlat2xy
    export latlon2yx_ja, lonlat2xy_ja
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using Printf

bin(x) = [xi > 0.5 for xi in x]

"""
    fix!(instance)

Verifies that the given unit commitment instance is valid and automatically fixes
some validation errors if possible, issuing a warning for each error found.
If a validation error cannot be automatically fixed, issues an exception.

Returns the number of validation errors found.
"""
function fix!(instance::UnitCommitmentInstance)::Int
    n_errors = 0

    for g in instance.units

        # Startup costs and delays must be increasing
        for s in 2:length(g.startup_categories)
            if g.startup_categories[s].delay <= g.startup_categories[s-1].delay
                prev_value = g.startup_categories[s].delay
                new_value = g.startup_categories[s-1].delay + 1
                @warn "Generator $(g.name) has non-increasing startup delays (category $s). " *
                      "Changing delay: $prev_value → $new_value"
                g.startup_categories[s].delay = new_value
                n_errors += 1
            end

            if g.startup_categories[s].cost < g.startup_categories[s-1].cost
                prev_value = g.startup_categories[s].cost
                new_value = g.startup_categories[s-1].cost
                @warn "Generator $(g.name) has decreasing startup cost (category $s). " *
                      "Changing cost: $prev_value → $new_value"
                g.startup_categories[s].cost = new_value
                n_errors += 1
            end

        end

        for t in 1:instance.time
            # Production cost curve should be convex
            for k in 2:length(g.cost_segments)
                cost = g.cost_segments[k].cost[t]
                min_cost = g.cost_segments[k-1].cost[t]
                if cost < min_cost - 1e-5
                    @warn "Generator $(g.name) has non-convex production cost curve " *
                          "(segment $k, time $t). Changing cost: $cost → $min_cost"
                    g.cost_segments[k].cost[t] = min_cost
                    n_errors += 1
                end
            end

            # Startup limit must be greater than min_power
            if g.startup_limit < g.min_power[t]
                new_limit = g.min_power[t]
                prev_limit = g.startup_limit
                @warn "Generator $(g.name) has startup limit lower than minimum power. " *
                      "Changing startup limit: $prev_limit → $new_limit"
                g.startup_limit = new_limit
                n_errors += 1
            end
        end
    end

    return n_errors
end

function validate(instance_filename::String, solution_filename::String)
    instance = UnitCommitment.read(instance_filename)
    solution = JSON.parse(open(solution_filename))
    return validate(instance, solution)
end

"""
    validate(instance, solution)::Bool

Verifies that the given solution is feasible for the problem. If feasible,
silently returns true. In infeasible, returns false and prints the validation
errors to the screen.

This function is implemented independently from the optimization model in `model.jl`, and
therefore can be used to verify that the model is indeed producing valid solutions. It
can also be used to verify the solutions produced by other optimization packages.
"""
function validate(instance::UnitCommitmentInstance,
                  solution::Union{Dict,OrderedDict};
                 )::Bool
    err_count = 0
    err_count += validate_units(instance, solution)
    err_count += validate_reserve_and_demand(instance, solution)

    if err_count > 0
        @error "Found $err_count validation errors"
        return false
    end

    return true
end

function validate_units(instance, solution; tol=0.01)
    err_count = 0

    for unit in instance.units
        production = solution["Production (MW)"][unit.name]
        reserve = solution["Reserve (MW)"][unit.name]
        actual_production_cost = solution["Production cost (\$)"][unit.name]
        actual_startup_cost = solution["Startup cost (\$)"][unit.name]
        is_on = bin(solution["Is on"][unit.name])
        switch_off = bin(solution["Switch off"][unit.name]) # some formulations may not use this

        for t in 1:instance.time
            # Auxiliary variables
            if t == 1
                is_starting_up = (unit.initial_status < 0) && is_on[t]
                is_shutting_down = (unit.initial_status > 0) && !is_on[t]
                ramp_up = max(0, production[t] + reserve[t] - unit.initial_power)
                ramp_down = max(0, unit.initial_power - production[t])
            else
                is_starting_up = !is_on[t-1] && is_on[t]
                is_shutting_down = is_on[t-1] && !is_on[t]
                ramp_up = max(0, production[t] + reserve[t] - production[t-1])
                #ramp_down = max(0, production[t-1] - production[t])
                ramp_down = max(0, production[t-1] + reserve[t-1] - production[t])
            end

            # Compute production costs
            production_cost, startup_cost = 0, 0
            if is_on[t]
                production_cost += unit.min_power_cost[t]
                residual = max(0, production[t] - unit.min_power[t])
                for s in unit.cost_segments
                    cleared = min(residual, s.mw[t])
                    production_cost += cleared * s.cost[t]
                    residual = max(0, residual - s.mw[t])
                end
            end

            # Production should be non-negative
            if production[t] < -tol
                @error @sprintf("Unit %s produces negative amount of power at time %d (%.2f)",
                                unit.name, t, production[t])
                err_count += 1
            end

            # Verify must-run
            if !is_on[t] && unit.must_run[t]
                @error @sprintf("Must-run unit %s is offline at time %d",
                                unit.name, t)
                err_count += 1
            end

            # Verify reserve eligibility
            if !unit.provides_spinning_reserves[t] && reserve[t] > tol
                @error @sprintf("Unit %s is not eligible to provide spinning reserves at time %d",
                                unit.name, t)
                err_count += 1
            end

            # If unit is on, must produce at least its minimum power
            if is_on[t] && (production[t] < unit.min_power[t] - tol)
                @error @sprintf("Unit %s produces below its minimum limit at time %d (%.2f < %.2f)",
                                unit.name, t, production[t], unit.min_power[t])
                err_count += 1
            end

            # If unit is on, must produce at most its maximum power
            if is_on[t] && (production[t] + reserve[t] > unit.max_power[t] + tol)
                @error @sprintf("Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)",
                                unit.name, t, production[t], reserve[t], unit.max_power[t])
                err_count += 1
            end

            # If unit is off, must produce zero
            if !is_on[t] && production[t] + reserve[t] > tol
                @error @sprintf("Unit %s produces power at time %d while off",
                                unit.name, t)
                err_count += 1
            end

            # Startup limit
            if is_starting_up && (ramp_up > unit.startup_limit + tol)
                @error @sprintf("Unit %s exceeds startup limit at time %d (%.2f > %.2f)",
                                unit.name, t, ramp_up, unit.startup_limit)
                err_count += 1
            end

            # Shutdown limit
            if is_shutting_down && (ramp_down > unit.shutdown_limit + tol)
              @error @sprintf("Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)\n\tproduction[t-1] = %.2f\n\treserve[t-1] = %.2f\n\tproduction[t] = %.2f\n\treserve[t] = %.2f\n\tis_on[t-1] = %d\n\tis_on[t] = %d",
                                unit.name, t, ramp_down, unit.shutdown_limit,
                                (t == 1 ? unit.initial_power : production[t-1]), production[t],
                                (t == 1 ? 0. : reserve[t-1]), reserve[t],
                                (t == 1 ? unit.initial_status != nothing && unit.initial_status > 0 : is_on[t-1]), is_on[t]
                               )
                err_count += 1
            end

            # Ramp-up limit
            if !is_starting_up && !is_shutting_down && (ramp_up > unit.ramp_up_limit + tol)
                @error @sprintf("Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)",
                                unit.name, t, ramp_up, unit.ramp_up_limit)
                err_count += 1
            end

            # Ramp-down limit
            if !is_starting_up && !is_shutting_down && (ramp_down > unit.ramp_down_limit + tol)
                @error @sprintf("Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)\n\tproduction[t-1] = %.2f\n\treserve[t-1] = %.2f\n\tproduction[t] = %.2f\n\treserve[t] = %.2f\n\tis_on[t-1] = %d\n\tis_on[t] = %d",
                                unit.name, t, ramp_down, unit.ramp_down_limit,
                                (t == 1 ? unit.initial_power : production[t-1]), production[t],
                                (t == 1 ? 0. : reserve[t-1]), reserve[t],
                                (t == 1 ? unit.initial_status != nothing && unit.initial_status > 0 : is_on[t-1]), is_on[t]
                               )
                err_count += 1
            end

            # Verify startup costs & minimum downtime
            if is_starting_up

                # Calculate how much time the unit has been offline
                time_down = 0
                for k in 1:(t-1)
                  if !is_on[t - k]
                        time_down += 1
                    else
                        break
                    end
                end
                if t == time_down + 1 && !switch_off[1]
                    # If unit has always been off, then the correct startup cost depends on how long was it off before t = 1
                    # Absent known initial conditions, we assume it was off for the minimum downtime
                    # TODO: verify the formulations are making the same assumption...
                    initial_down = unit.min_downtime
                    if unit.initial_status < 0
                        initial_down = -unit.initial_status
                    end
                    time_down += initial_down
                end

                # Calculate startup costs
                for c in unit.startup_categories
                    if time_down >= c.delay
                        startup_cost = c.cost
                    end
                end

                # Check minimum downtime
                if time_down < unit.min_downtime
                    @error @sprintf("Unit %s violates minimum downtime at time %d",
                                    unit.name, t)
                    err_count += 1
                end
            end

            # Verify minimum uptime
            if is_shutting_down

                # Calculate how much time the unit has been online
                time_up = 0
                for k in 1:(t-1)
                    if is_on[t - k]
                        time_up += 1
                    else
                        break
                    end
                end
                if t == time_up + 1
                    initial_up = unit.min_uptime
                    if unit.initial_status > 0
                        initial_up = unit.initial_status
                    end
                    time_up += initial_up
                end

                if (t == time_up + 1) && (unit.initial_status > 0)
                    time_up += unit.initial_status
                end

                # Check minimum uptime
                if time_up < unit.min_uptime
                    @error @sprintf("Unit %s violates minimum uptime at time %d",
                                    unit.name, t)
                    err_count += 1
                end
            end

            # Verify production costs
            if abs(actual_production_cost[t] - production_cost) > 1.00
                @error @sprintf("Unit %s has unexpected production cost at time %d (%.2f should be %.2f)",
                                unit.name, t, actual_production_cost[t], production_cost)
                err_count += 1
            end

            # Verify startup costs
            if abs(actual_startup_cost[t] - startup_cost) > 1.00
                @error @sprintf("Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)",
                                unit.name, t, actual_startup_cost[t], startup_cost)
                err_count += 1
            end

        end
    end

    return err_count
end

function validate_reserve_and_demand(instance, solution, tol=0.01)
    err_count = 0
    for t in 1:instance.time
        load_curtail = 0
        fixed_load = sum(b.load[t] for b in instance.buses)
        production = sum(solution["Production (MW)"][g.name][t]
                         for g in instance.units)
        if "Load curtail (MW)" in keys(solution)
            load_curtail = sum(solution["Load curtail (MW)"][b.name][t]
                               for b in instance.buses)
        end
        balance = fixed_load - load_curtail - production

        # Verify that production equals demand
        if abs(balance) > tol
            @error @sprintf("Non-zero power balance at time %d (%.2f - %.2f - %.2f != 0)",
                            t, fixed_load, load_curtail, production)
            err_count += 1
        end

        # Verify spinning reserves
        reserve = sum(solution["Reserve (MW)"][g.name][t] for g in instance.units)
        if reserve < instance.reserves.spinning[t] - tol
            @error @sprintf("Insufficient spinning reserves at time %d (%.2f should be %.2f)",
                            t, reserve, instance.reserves.spinning[t])
            err_count += 1
        end
    end

    return err_count
end


text:
function uncertainty_set_with_index(i, data, ϵ)
    set = uncertainty_set(data[i, :], ϵ)
    return collect(zip(fill(i, length(set)), set))
end

function combine(a::NTuple{K, Tuple{CartesianIndex{M}, Vector{Float64}}}, s::NTuple{N, Int})::Array{Float64, N} where {M, N, K}
    @assert M == N - 1
    p = zeros(s...)
    for (i, v) in a
        p[i, :] = v
    end
    return p
end

function uncertainty_sets(data::Array{Float64, N}, ϵ::Float64)::Vector{Array{Float64, N}} where N
    # TODO: ϵs::Vector{Float64}, size: |S_j|, j::Node
    s = size(data)
    sets = [uncertainty_set_with_index(i, data, ϵ) for i in CartesianIndices(s[1:end-1])];
    prods = vec(map(a -> combine(a, s), product(sets...)));
    return prods
end

function worst_case(model, objs)
    @variable(model, x)
    for obj in objs
        @constraint(model, x ≤ obj)
    end
    return x
end

max_stars_repo_name: jaantollander/RobustDecisionProgramming

max_issues_repo_name: jaantollander/RobustDecisionProgramming

max_forks_repo_name: jaantollander/RobustDecisionProgramming


text:
macro rand(exp)
    rand_macro(exp)
end

function rand_macro(ex)
    whereparams = []
    ex isa Expr && ex.head ∈ (:(=), :function, :->) ||
        throw(ArgumentError("@rand requires an expression defining `rand`"))
    sig = ex.args[1]
    body = ex.args[2]

    if sig.head == :where
        append!(whereparams, sig.args[2:end])
        sig = sig.args[1]
    end

    if ex.head == :function && sig.head == :tuple # anonymous function
        sig = Expr(:call, :rand, sig.args...)
    end
    if ex.head == :->
        # TODO: check that only one argument is passed
        sig = Expr(:call, :rand, sig)
    end

    sig.head == :call &&
        sig.args[1] == :rand || throw(ArgumentError(
            "@rand requires a function expression defining `rand`"))

    argname = sig.args[2].args[1] # x
    namefull = sig.args[2] # x::X
    @assert namefull.head == :(::) # TODO: throw exception

    # sub-samplers; second argument true forces Val(Inf) for the sampler
    sps = Pair{<:Any,Bool}[]

    rng = gensym()
    body = samplerize!(sps, body, argname, rng)
    istrivial = isempty(sps)

    exsig = Expr(:call,
                 :(Random.rand),
                 :($(esc(rng))::AbstractRNG),
                 esc(as_sampler(namefull, istrivial)))

    if !isempty(whereparams)
        exsig = Expr(:where, exsig, map(esc, whereparams)...)
    end

    ex = Expr(:function, exsig, esc(body))

    sp = if istrivial
        # we explicitly define Sampler even in the trivial case to handle
        # redefinitions, where the old rand/sampler pair (for SamplerSimple)
        # is overwritten by a new one (for SamplerTrivial)
        quote
            Random.Sampler(::Type{RNG}, n::Repetition) where {RNG<:AbstractRNG} =
                SamplerTrivial($(esc(argname)))
        end
    else
        quote
            Random.Sampler(::Type{RNG}, n::Repetition) where {RNG<:AbstractRNG} =
                SamplerSimple($(esc(argname)), tuple(SP))
        end
    end
    @assert sp.args[2].args[1].head == :where
    append!(sp.args[2].args[1].args, map(esc, whereparams))

    # insert x::X in the argument list, between RNG and n::Repetition
    insert!(sp.args[2].args[1].args[1].args, 3, esc(namefull))

    # insert inner samplers
    if !istrivial
        SP = [Expr(:call, :Sampler, :RNG, esc(x), many ? Val{Inf}() : :n) for (x, many) in sps]
        @assert :SP == pop!(sp.args[2].args[2].args[2].args[3].args)
        append!(sp.args[2].args[2].args[2].args[3].args, SP)
    end

    quote
        $ex
        $(sp.args[2]) # unwrap the quote/block around the definition
    end
end

function as_sampler(ex, istrivial)
    t = istrivial ? :(RandomExtensions.SamplerTrivial) : :(RandomExtensions.SamplerSimple)
    Expr(:(::),
         ex.args[1],
         Expr(:curly, t,
              Expr(:(<:), ex.args[2])))
end

function samplerize!(sps, ex, name, rng)
    if ex == name
        # not within a rand() call
        return Expr(:ref, name) # name -> name[]
    end
    ex isa Expr || return ex
    if ex.head == :call && ex.args[1] == :rand
        # we assume that if rand has more than one arg, we want
        # a Val(Inf) sampler (e.g. rand(1:9, 2, 3)
        push!(sps, ex.args[2] => length(ex.args) > 2)
        i = length(sps)
        Expr(:call, :rand, rng, :($name.data[$i]), ex.args[3:end]...)
    else
        Expr(ex.head, map(e -> samplerize!(sps, e, name, rng), ex.args)...)
    end
end


text:
"""
$(TYPEDEF)

Bathe's integration scheme with sub-steps of equal size.

### Fields

- `Δt` -- step-size
- `α`  -- parameter α of the method
- `δ`  -- parameter δ of the method

### References

See [[BAT07]](@ref).
"""
struct Bathe{N} <: AbstractSolver
    Δt::N
end

Bathe(; Δt::N) where {N} = Bathe(Δt)

step_size(alg::Bathe) = alg.Δt

function _init(alg::Bathe, M, C, K)
    Δt = alg.Δt

    # compute integration constants
    ξ = 1/Δt
    ξ² = ξ*ξ
    a₀ = 16*ξ²
    a₁ = 4*ξ
    a₂ = 9*ξ²
    a₃ = 3*ξ
    a₄ = 2*a₁
    a₅ = 12*ξ²
    a₆ = -3*ξ²
    a₇ = -ξ

    K̂₁ = K + a₀ * M + a₁ * C
    K̂₂ = K + a₂ * M + a₃ * C
    return a₀, a₁, a₂, a₃, a₄, a₅, a₆, a₇, K̂₁, K̂₂
end

# a comment about notation: variables with superscript `+` correspond to
# evaluation at intermediate times, i.e. t + Δt/2
function _solve(alg::Bathe{N},
                ivp::InitialValueProblem{ST, XT},
                NSTEPS::Int; kwargs...) where {N, VT, ST, XT<:Tuple{VT, VT}}

    sys = system(ivp)
    (U₀, U₀′) = initial_state(ivp)

    IMAX = NSTEPS + 1
    M, C, K, R = _unwrap(sys, IMAX)

    U₀′′ = M \ (R[1] - C * U₀′ - K * U₀)
    a₀, a₁, a₂, a₃, a₄, a₅, a₆, a₇, K̂₁, K̂₂ = _init(alg, M, C, K)
    K̂₁⁻¹ = factorize(K̂₁)
    K̂₂⁻¹ = factorize(K̂₂)

    # initialize displacements, velocities and accelerations
    U = Vector{VT}(undef, IMAX)
    U′ = Vector{VT}(undef, IMAX)
    U′′ = Vector{VT}(undef, IMAX)
    U[1] = U₀
    U′[1] = U₀′
    U′′[1] = U₀′′

    @inbounds for i in 1:NSTEPS
        # ---------------------------------------------
        # First sub-step (Newmark trapezoidal rule)
        # ---------------------------------------------

        # calculate effective loads at time t + Δt/2
        mᵢ = M * (a₀ * U[i] + a₄ * U′[i] + U′′[i])
        cᵢ = C * (a₁ * U[i] + U′[i])
        Rᵢ⁺ = (R[i] + R[i+1])/2
        R̂ᵢ⁺ = Rᵢ⁺ + mᵢ + cᵢ

        # solve for displacements at time t + Δt/2
        Uᵢ⁺ = K̂₁⁻¹ \ R̂ᵢ⁺

        # calculate velocities at time t + Δt/2
        U′ᵢ⁺ = a₁ * (Uᵢ⁺ - U[i]) - U′[i]

        # ---------------------------------------------
        # Second sub-step (3-point Euler backward method)
        # ---------------------------------------------

        # calculate effective loads at time t + Δt
        mᵢ = M * (a₅ * Uᵢ⁺ + a₆ * U[i] + a₁ * U′ᵢ⁺ + a₇ * U′[i])
        cᵢ = C * (a₁ * Uᵢ⁺ + a₇ * U[i])
        R̂ᵢ₊₁ = R[i+1] + mᵢ + cᵢ

        # solve for displacements at time t + Δt
        U[i+1] = K̂₂⁻¹ \ R̂ᵢ₊₁

        # calculate velocities and accelerations at time t + Δt
        U′[i+1] = -a₇ * U[i] - a₁ * Uᵢ⁺ + a₃ * U[i+1]
        U′′[i+1] = -a₇ * U′[i] - a₁ * U′ᵢ⁺ + a₃ * U′[i+1]
    end

    return _build_solution(alg, U, U′, U′′, NSTEPS)
end

function _build_solution(alg::Bathe{N}, U, U′, U′′, NSTEPS) where {N}
    t = range(zero(N), step=alg.Δt, length=(NSTEPS+1))
    return Solution(alg, U, U′, U′′, t)
end


text:
using GeometricIntegrators.Integrators
using GeometricIntegrators.Integrators.SPARK
using GeometricIntegrators.Integrators.VPRK
using GeometricIntegrators.Tableaus
using GeometricIntegrators.Utils
using Test

include("coefficients_tests.jl")
include("vprk_tableaus_tests.jl")
include("spark_tableaus_tests.jl")
include("splitting_tableaus_tests.jl")


text:
using Contour, Test

include("verify_vertices.jl")
include("interface.jl")

@test length(detect_ambiguities(Base, Contour)) <= 4 # should be zero but there are a few stragglers in the imports

#issue 59
@inferred collect(())


text:
module execute2d
using LPDM
using POMDPs, Parameters, StaticArrays, D3Trees, Distributions, SparseArrays
using Combinatorics
using StatsFuns
using Printf
using Dates

include("LightDarkLPDM.jl")
using LightDarkPOMDPs

include("LPDMBounds2d.jl")

struct LPDMTest
    solver_mode::Symbol
    action_mode::Symbol
    obs_mode::Symbol
    reward_mode::Symbol
end

struct LPDMScenario
    s0::LD2State
end

# reward function options - :quadratic or :fixed
function batch_execute(;n::Int64=1, debug::Int64=1, reward_mode=:quadratic)
    test=Array{LPDMTest}(undef,0)

    push!(test, LPDMTest(:lpdm, :standard, :discrete, reward_mode))
    push!(test, LPDMTest(:lpdm, :standard, :continuous, reward_mode))
    # push!(test, LPDMTest(:lpdm_bv, :bv, reward_func)) # blind value
    # push!(test, LPDMTest(:lpdm, :adapt, reward_func)) # simulated annealing

    scen=Array{LPDMScenario}(undef,0)
    # push!(scen, LPDMScenario(LD2State(-2*π, π)))
    push!(scen, LPDMScenario(LD2State(-π, π)))
    # push!(scen, LPDMScenario(LD2State(π/2, π/2)))
    # push!(scen, LPDMScenario(LD2State(π, -π)))
    # push!(scen, LPDMScenario(LD2State(2*π,2*π)))

    # Dummy execution, just to make sure all the code is compiled and loaded,
    # to improve uniformity of subsequent executions.
    # execute(solv_mode = :despot, action_space_type = :small, n_sims = 1, s0 = LD2State(π,π), output = 0)
    # execute(solv_mode = :lpdm, action_space_type = :adaptive, n_sims = 1, s0 = LD2State(π,π), output = 0)

    # General solver parameters
    steps::Int64                = -1
    time_per_move::Float64      = 1.0
    search_depth::Int64         = 30
    n_particles::Int64          = 50
    max_trials::Int64           = -1

    f = open("results_" * Dates.format(now(),"yyyy-mm-dd_HH_MM") * ".txt", "w")

    Printf.@printf(f,"GENERAL SOLVER PARAMETERS\n")
    Printf.@printf(f,"\treward function:\t\t\t%s\n", reward_func)
    Printf.@printf(f,"\tsteps:\t\t\t%d\n", steps)
    Printf.@printf(f,"\ttime per move:\t\t\t%f\n", time_per_move)
    Printf.@printf(f,"\tsearch depth:\t\t\t%d\n", search_depth)
    Printf.@printf(f,"\tN particles:\t\t\t%d\n", n_particles)
    Printf.@printf(f,"\tmax trials:\t\t\t%d\n\n", max_trials)

    for i in 1:length(scen)
        # write(f,"SCENARIO $i, s0 = $(scen[i].s0)\n\n")
        if debug >= 0
            println("")
            println("SCENARIO $i, s0 = $(scen[i].s0)")
            println("------------------------")
        end

        Printf.@printf(f,"SCENARIO %d, s0 = (%f,%f)\n", i, scen[i].s0[1], scen[i].s0[2])
        Printf.@printf(f,"==================================================================\n")
        # Printf.@printf(f,"mode\t\tact. space\t\tsteps(std)\t\treward(std)\n")
        Printf.@printf(f,"SOLVER\t\tACT. MODE\t\tOBS. MODE\t\tSTEPS (STD)\t\t\tREWARD (STD)\n")
        Printf.@printf(f,"=====================================================================\n")
        for t in test
            if debug >= 0
                println("mode: $(t.solver_mode), actions: $(t.action_mode), observations: $(t.obs_mode), rewards: $(t.reward_mode)")
            end
            steps_avg, steps_std, reward_avg, reward_std =
                        execute(solv_mode         = t.mode,
                                action_space_type = t.action_space,
                                n_sims            = n,
                                steps             = steps,
                                time_per_move     = time_per_move,
                                search_depth      = search_depth,
                                n_particles       = n_particles,
                                max_trials        = max_trials,
                                s0                = scen[i].s0,
                                output            = debug,
                                reward_func       = reward_func)

            Printf.@printf(f,"%s\t\t%s\t\t\t%05.2f (%06.2f)\t\t%06.2f (%06.2f)\n",
                            string(t.mode), string(t.action_space), steps_avg, steps_std, reward_avg, reward_std)

            debug >=0 && println("STATS: solver=$(t.mode), actions=$(t.action_space), steps = $(steps_avg) ($steps_std), reward = $(reward_avg) ($reward_std)")
        end
        Printf.@printf(f,"==================================================================\n")
        Printf.@printf(f,"%d tests per scenario\n\n", n)
    end
    close(f)
end

function execute(;vis::Vector{Int64}        = Int64[],
                solv_mode::Symbol           = :despot,
                action_space_type::Symbol   = :small,
                reward_func::Symbol         = :quadratic,
                n_sims::Int64               = 1,
                steps::Int64                = -1,
                time_per_move::Float64      = 5.0,
                search_depth::Int64         = 50,
                n_particles::Int64          = 50,
                max_trials::Int64           = -1,
                # s0::LD2State                = LD2State(0.5*π, 2*π),
                s0::LD2State                = LD2State(π, -π),
                output::Int64               = 1
                )

    if solv_mode == :despot
        p = LightDark2DDespot(action_space_type, reward_func = reward_func)
    elseif (solv_mode == :lpdm || solv_mode == :lpdm_bv)
        p = LightDark2DLpdm(action_space_type, reward_func = reward_func)
    end

    sim_rewards = Vector{Float64}(undef,n_sims)
    sim_steps   = Vector{Int64}(undef,n_sims)

    for sim in 1:n_sims

        world_rng = RNGVector(1, UInt32(sim))
        LPDM.set!(world_rng, 1)

        # NOTE: restrict initial state to positive numbers only, for now
        # s::LD2State                  = LD2State(π);    # initial state
        step_rewards::Array{Float64}     = Vector{Float64}(undef,0)

        # NOTE: Original version
        # solver = LPDM.LPDMSolver{LD2State, LD2Action, LD2Obs, LDBounds2d{LD2State, LD2Action, LD2Obs}, RNGVector}(
        #                                                                     # rng = sim_rng,
        #                                                                     debug = output,
        #                                                                     time_per_move = 1.0,  #sec
        #                                                                     # time_per_move = 1.0,  #sec
        #                                                                     sim_len = -1,
        #                                                                     search_depth = 50,
        #                                                                     n_particles = 20,
        #                                                                     # seed = UInt32(5),
        #                                                                     seed = UInt32(2*sim+1),
        #                                                                     # max_trials = 10)
        #                                                                     max_trials = -1,
        #                                                                     mode = solv_mode)

        # NOTE: Test version
        solver = LPDM.LPDMSolver{LD2State, LD2Action, LD2Obs, LDBounds2d{LD2State, LD2Action, LD2Obs}, RNGVector}(
                                                                            # rng = sim_rng,
                                                                            debug = output,
                                                                            time_per_move = time_per_move,  #sec
                                                                            # time_per_move = 1.0,  #sec
                                                                            sim_len = steps,
                                                                            #sim_len = 20,
                                                                            search_depth = search_depth,
                                                                            n_particles = n_particles,
                                                                            # seed = UInt32(2),
                                                                            seed = UInt32(2*sim+1),
                                                                            # max_trials = 1000)
                                                                            max_trials = max_trials,
                                                                            mode = solv_mode)

        # #DEBUG: remove ###############################
        # # states=[Vec2(3.6971,7.3388), Vec2(3.5,-1.2), Vec2(7.1,7.3), Vec2(5.35,0.0)]
        # states=[Vec2(3.14,3.14)]
        # for st in states
        #     b = LDBounds2d{LD2State, LD2Action, LD2Obs}(p)
        #     lb, ub = LPDM.bounds(b,p,[LPDMParticle{Vec2}(st,1,1.0)],solver.config)
        #     println("s=$st, lb = $lb, ub = $ub, lba=$(best_lb_action(b)), uba=$(best_ub_action(b))")
        # end
        # error("stop for now")
        # ##############################################
    #---------------------------------------------------------------------------------
        # Belief
        bu = LPDMBeliefUpdater(p,
                               n_particles = solver.config.n_particles,
                               seed = UInt32(3*sim+1));  # initialize belief updater
        initial_states = state_distribution(p, s0, solver.config, world_rng)     # create initial  distribution
        # println(initial_states)
        current_belief = LPDM.create_belief(bu)                       # allocate an updated belief object
        LPDM.initialize_belief(bu, initial_states, current_belief)    # initialize belief
        # show(current_belief)
        updated_belief = LPDM.create_belief(bu)
        # println("$(length(actions(p))) actions: $(actions(p))")
        # error("")
    #---------------------------------------------------------------------------------

        policy::LPDMPolicy = POMDPs.solve(solver, p)
        s = s0
        step::Int64 = 0
        r::Float64 = 0.0

        # output >= 1 && println("=============== SIMULATION # $sim ================")
        if output >= 0
            println("")
            println("*** SIM $sim (mode: $solv_mode, action space: $action_space_type, s0: $s0)***")
            println("")
        end

        val, run_time, bytes, gctime, memallocs =
        @timed while !isterminal(p, s) && (solver.config.sim_len == -1 || step < solver.config.sim_len)
            step += 1

            if output >= 1
                println("")
                println("=============== Step $step ================")
                show(current_belief)
            end

            a = POMDPs.action(policy, current_belief)
            if output >= 1
                println("s: $s")
                println("a: $a")
            end

            s, o, r = POMDPs.generate_sor(p, s, a, world_rng)
            push!(step_rewards, r)
            if output >= 1
                println("s': $s")
                println("o': $o")
                println("r: $r")
                println("=======================================")
            end

            # update belief
            POMDPs.update(bu, current_belief, a, o, updated_belief)
            current_belief = deepcopy(updated_belief)
            # show(updated_belief) #NOTE: don't show for now

            if LPDM.isterminal(p, current_belief.particles)
                if output >= 1
                    println("Terminal belief. Execution completed.")
                    show(current_belief)
                end
                break
            end

            if step ∈ vis
                t = LPDM.d3tree(solver,
                                detect_repeat=false,
                                title="Step $step",
                                init_expand=10)
                # # show(t)
                inchrome(t)
                # blink(t)
            end
            if output >= 1
                solv_mode == :lpdm && println("root actions: $(solver.root.action_space)")
            end
        end

        sim_steps[sim] = step
        sim_rewards[sim] = sum(step_rewards)
        output >= 0 && println("steps=$(sim_steps[sim]), reward=$(sim_rewards[sim])")
    end

    return mean(sim_steps), std(sim_steps), mean(sim_rewards), std(sim_rewards)

end

end #module

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
using EmpiricalBayes
using Test
using Distributions
using StatsBase

true_dist = MixtureModel([ Normal(-0.3,.5), Normal(1.05,.5)])
marginal_grid = collect(range(-6, stop=6, length=1001));

marginal_h = marginal_grid[2] - marginal_grid[1] # should really replace by range

prior_grid = collect(range(-3, stop=3, length=121));
d_true = NormalConvolutionProblem(true_dist, marginal_grid);

ds = MixingNormalConvolutionProblem(Normal, 0.2, prior_grid, marginal_grid);

m_train= 10_000

Xs_train = rand(d_true,m_train)
Xs_test = rand(d_true,m_train)

target = PosteriorTarget(PosteriorMeanNumerator(1.0))

M_bd = marginal_h/sqrt(2*pi)
f_const = BinnedMarginalDensity(M_bd, marginal_grid, marginal_h)
# Train
# Estimate the numerator
M_max_num = MinimaxCalibrator(ds, f_const, m_train, target.num; C=Inf)
M_max_denom = MinimaxCalibrator(ds, f_const, m_train, target.denom; C=Inf)

ceb_ci_1 =  CEB_ci(Xs_train, Xs_test,
            ds,
            target,
            M_max_num, #
            M_max_denom)

# Check if CB intervals work

CB_num = ComteButucea(target.num, m_train, marginal_grid)
CB_denom = ComteButucea(target.denom, m_train, marginal_grid)

ceb_ci_1_cb =  CEB_ci(Xs_train, Xs_test,
            ds,
            target,
            CB_num, #
            CB_denom)

# Both approaches perform remarkably similar in this case
posterior_stats(d_true, target)[3] #true target
confint(ceb_ci_1, target)
confint(ceb_ci_1_cb, target)

ceb_ci_1_cb.est_target, ceb_ci_1.est_target

brad_g =  fit(BradDeconvolveR, Xs_train, prior_grid=range(-3.6, stop=3.6, length=361);
    deltaAt = :nothing, pDegree=10)

ceb_ci_brad = CEB_ci(Xs_train, Xs_test,ds, target, brad_g)

posterior_stats(d_true, target)[3]

confint(ceb_ci_brad, target)
ceb_ci_brad.est_target
ceb_ci_brad.calibrated_target


text:
using PrimeSieve
using Base.Test

include("msieve_test.jl")
include("gaussian_test.jl")
include("tuples_test.jl")
include("primecount_test.jl")
include("primesieve_test.jl")
include("test_readme.jl")  # this is slow too.
include("slow_test.jl")


text:
using Test, SummationByPartsOperators

for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
    xmin = T(-1)
    xmax = T(1)
    N = 2^6
    u0func = sinpi
    tspan = (zero(T), one(T))
    afunc = one

    # Legendre
    let D = legendre_derivative_operator(xmin, xmax, N)
        Di = nothing
        semidisc = VariableLinearAdvectionNonperiodicSemidiscretisation(D, Di, afunc, split_form, zero, zero)
        ode = semidiscretise(u0func, semidisc, tspan)
        du = similar(ode.u0)
        semidisc(du, ode.u0, nothing, first(tspan))
    end

    # SBP FD
    for acc_order in 2:2:8
        D = derivative_operator(MattssonSvärdNordström2004(), 1, acc_order, xmin, xmax, N)
        Di = dissipation_operator(D)
        semidisc = VariableLinearAdvectionNonperiodicSemidiscretisation(D, Di, afunc, split_form, zero, zero)
        ode = semidiscretise(u0func, semidisc, tspan)
        du = similar(ode.u0)
        semidisc(du, ode.u0, nothing, first(tspan))
    end
end