text:
using Pkg
Pkg.activate(pwd())

using titanic

traindata, trainlabels = loadtrain("../data/train.csv"; normalize=false)
testdata = loadtest("../data/test.csv"; normalize=false)

H = boost(traindata, trainlabels)

prediction = classify(testdata, H)

savepredictions("../results/predictions-adaboost.csv", prediction)

ext: jl


text:
function EMPolya(alpha0::AbstractArray,delta::AbstractArray,alpha::AbstractArray,Y::AbstractArray,TS::AbstractArray,p::AbstractArray,L::Int64,I::AbstractArray,J::Int64,K::Int64,EMNum::Int64,MMNum::Int64,BaseID::Int64=0,stopc::Float64=1e-4)
    mu=Array{Array{Float64}}(undef,L)
    MTS=zeros(Int64,L)
    MY=zeros(Int64,J,L)
    dLike=0
    dLikeNew=0
    MLike=zeros(Float64,K)
    lp=log.(p)
    prog = ProgressThresh(stopc, "Criterion:")

    #Average the initial delta
    if L>1
        delta[:,:,1:end.!=BaseID].=mean(delta[:,:,1:end.!=BaseID],dims=3)
        alpha=alpha0.+delta
    end

        @inbounds    for l in 1:L
            mu[l]=zeros(Float64,K,I[l])
            MTS[l]=maximum(TS[l])
            @inbounds    for i in 1:I[l]
                @inbounds for k in 1:K
                    MLike[k]=lp[k,l]+logpdf(DirichletMultinomial(TS[l][i],alpha[:,k,l]),view(Y[l],:,i))
                end
                m=maximum(MLike)
                @inbounds for k in 1:K
                    mu[l][k,i]=MLike[k]-m-log(sum(exp.(MLike[:].-m)))
                end
                dLike+=m+log(sum(exp.(MLike[:].-m)))

        end
        lp[:,l]=maximum(mu[l],dims=2).+log.(sum(exp.(mu[l].-maximum(mu[l],dims=2)),dims=2)).-log(I[l])
        @inbounds for j in 1:J
            MY[j,l]=maximum(Y[l][j,:])
        end
    end
    Newmu=copy(mu)
    Newlp=copy(lp)

    Flag=1
    @inbounds for t1 in 1:EMNum

        #For the first ten EM iterations, Recalculate and average the delta
        if 2<=t1<=10
            if L>1
                delta[:,:,1:end.!=BaseID].=mean(delta[:,:,1:end.!=BaseID],dims=3)
                alpha=alpha0.+delta
            end

            @inbounds    for l in 1:L
                @inbounds    for i in 1:I[l]
                    @inbounds for k in 1:K
                        MLike[k]=lp[k,l]+logpdf(DirichletMultinomial(TS[l][i],alpha[:,k,l]),view(Y[l],:,i))
                    end
                    m=maximum(MLike)
                    @inbounds for k in 1:K
                        Newmu[l][k,i]=MLike[k]-m-log(sum(exp.(MLike[:].-m)))
                    end
                    dLike+=m+log(sum(exp.(MLike[:].-m)))

            end

            end

        end

        #MM step
        alpha0, delta, alpha=MMPolya(alpha0,delta,alpha,Y,TS,MY,MTS,Newmu,L,I,J,K,MMNum)

        #Compute the updated likelihood
        dLikeNew=0
        @inbounds for l in 1:L
            @inbounds for i in 1:I[l]
                s=0
                @inbounds for k in 1:K
                    MLike[k]=lp[k,l]+logpdf(DirichletMultinomial(TS[l][i],alpha[:,k,l]),view(Y[l],:,i))
                end
                m=maximum(MLike)
                @inbounds for k in 1:K
                    Newmu[l][k,i]=MLike[k]-m-log(sum(exp.(MLike[:].-m)))
                end
                dLikeNew+=m+log(sum(exp.(MLike[:].-m)))
            end
            Newlp[:,l]=maximum(mu[l],dims=2).+log.(sum(exp.(mu[l].-maximum(mu[l],dims=2)),dims=2)).-log(I[l])
        end

        #Stop criteria
        ProgressMeter.update!(prog, abs((dLikeNew-dLike)/dLike))
        sleep(0.5)

        if dLikeNew>dLike
            if abs((dLikeNew-dLike)/dLike)> stopc
                dLike=dLikeNew
                mu=copy(Newmu)
                lp=copy(Newlp)
            else
                Flag=0
            end
        else
            Flag=0
            println("decreasing likelihood!")
            println(dLike)
            println(dLikeNew)
        end

        if(Flag==0)
            break
        end

    end

    #Determine cell label based on mu
    groups=mapslices(argmax, mu[1], dims=1)[:]
    if L>1
        @inbounds for l in 2:L
            groups=vcat(groups,mapslices(argmax, mu[l], dims=1)[:])
        end
    end

    return groups, dLikeNew, alpha0, delta, alpha, lp
end


text:
function solve()
    n = parse(Int, readline())
    while n % 2 == 0
        n ÷= 2
    end
    sq = floor(Int, sqrt(n))
    ret = 0
    for i=1:sq
        if n % i == 0
            ret += 2
        end
    end
    if sq * sq == n
        ret -= 1
    end
    2ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Documenter
using MERAKit

makedocs(sitename = "MERAKit.jl", modules = [MERAKit])
deploydocs(repo = "github.com/mhauru/MERAKit.jl.git")


text:
using Hyperscript, Markdown
using JSServe, Observables
using JSServe: Session, evaljs, linkjs, div, active_sessions
using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, Asset
using WGLMakie, AbstractPlotting

markdown_css = Asset(JSServe.dependency_path("markdown.css"))

function test_handler(session, req)
    button = Button("click")
    slider = Slider(1:100)

    signal = map(slider) do value
        sin.(LinRange(0, value, 1000))
    end

    scene = scatter(signal, markersize=(10.0,10.0), resolution=(500,200), limits=FRect(0, -1, 1000, 2))

    on(button) do val
        scene[end].color = rand(RGBf0)
    end
    dom = md"""
    # More MD

    [Github-flavored Markdown info page](http://github.github.com/github-flavored-markdown/)

    [![Build Status](https://travis-ci.com/SimonDanisch/JSServe.jl.svg?branch=master)](https://travis-ci.com/SimonDanisch/JSServe.jl)

    Lalala
    ======

    Alt-H2
    ------

    *italic* or **bold**

    Combined emphasis with **asterisks and _underscores_**.

    1. First ordered list item
    2. Another item
        * Unordered sub-list.
    1. Actual numbers don't matter, just that it's a number
        1. Ordered sub-list

    * Unordered list can use asterisks

    Inline `code` has `back-ticks around` it.
    ```julia
    test("haha")
    ```

    ---
    # JSServe

    [![Build Status](https://travis-ci.com/SimonDanisch/JSServe.jl.svg?branch=master)](https://travis-ci.com/SimonDanisch/JSServe.jl)
    [![Build Status](https://ci.appveyor.com/api/projects/status/github/SimonDanisch/JSServe.jl?svg=true)](https://ci.appveyor.com/project/SimonDanisch/JSServe-jl)
    [![Codecov](https://codecov.io/gh/SimonDanisch/JSServe.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/SimonDanisch/JSServe.jl)
    [![Build Status](https://travis-ci.com/SimonDanisch/JSServe.jl.svg?branch=master)](https://travis-ci.com/SimonDanisch/JSServe.jl)

    | Tables        | Are           | Cool  |
    | ------------- |:-------------:| -----:|
    | col 3 is      | right-aligned | $1600 |
    | col 2 is      | centered      |   $12 |
    | zebra stripes | are neat      |    $1 |

    > Blockquotes are very handy in email to emulate reply text.
    > This line is part of the same quote.

    # Plots:

    $(slider)

    $(button)

    ---

    $(scene)

    ---
    """
    return DOM.div(markdown_css, dom)
end

app = JSServe.Application(test_handler, "0.0.0.0", 8081)


text:
using ApproximateBayesianComputing
const ABC = ApproximateBayesianComputing
using Distributions
using Random

Random.seed!(1234)

# Set Prior for Population Parameters
theta_true = [0.0, 1.0]
param_prior = Distributions.MvNormal(theta_true,ones(length(theta_true)))

# Code to generate simulated data given array of model parameters
num_data_default = 100
gen_data_normal(theta::Array, n::Integer = num_data_default) = rand(Normal(theta[1],theta[2]),num_data_default)

# Function to adjust originally proposed model parameters, so that they will be valid
normalize_theta2_pos(theta::Array) =  theta[2] = abs(theta[2])

# Function to test if the proposed model parameter are valid
is_valid_theta2_pos(theta::Array) =  theta[2]>0.0 ? true : false

# Tell ABC what it needs to know for a simulation
abc_plan = abc_pmc_plan_type(gen_data_normal,ABC.calc_summary_stats_mean_var,ABC.calc_dist_max, param_prior; is_valid=is_valid_theta2_pos,num_max_attempt=10000);

# Generate "true/observed data" and summary statistics
data_true = abc_plan.gen_data(theta_true)
ss_true = abc_plan.calc_summary_stats(data_true)
#println("theta= ",theta_true," ss= ",ss_true, " d= ", 0.)

# Run ABC simulation
@time pop_out = run_abc(abc_plan,ss_true;verbose=true);

#=
# Optional plotting of results
using PyPlot

hist(pop_out.weights*length(pop_out.weights));
hist(pop_out.dist);

num_param = 2
num_grid_x = 100
num_grid_y = 100
limit = 1.0
x = collect(linspace(theta_true[1]-limit,theta_true[1]+limit,num_grid_x));
y = collect(linspace(theta_true[2]-limit,theta_true[2]+limit,num_grid_y));
z = zeros(Float64,(num_param,length(x),length(y)))
for i in 1:length(x), j in 1:length(y)
    z[1,i,j] = x[i]
    z[2,i,j] = y[j]
end
z = reshape(z,(num_param,length(x)*length(y)))
zz = [ ABC.pdf(ABC.GaussianMixtureModelCommonCovar(pop_out.theta,pop_out.weights,ABC.cov_weighted(pop_out.theta',pop_out.weights)),vec(z[:,i])) for i in 1:size(z,2) ]
zz = reshape(zz ,(length(x),length(y)));
levels = [exp(-0.5*i^2)/sqrt(2pi^num_param) for i in 0:5];
PyPlot.contour(x,y,zz',levels);
plot(pop_out.theta[1,:],pop_out.theta[2,:],".");
=#


text:
1 + 2

x = 10
x + 1
x = 1 + 1
x = "Hello, World!"

pi

typeof(1)

Sys.WORD_SIZE

Int

UInt

0x1

0b10

0o010

(typemin(Int32), typemax(Int32))

1.0

5 % 2

x = 3

2x

[1,2,3] .+ 3

'x'

typeof('x')

greet = "Hello"
whom = "world"
"$greet, $whom.\n"

max_stars_repo_name: kpvla/til

max_issues_repo_name: kpvla/til

max_forks_repo_name: kpvla/til


text:
using Test, YaoBlocks, YaoArrayRegister
import YaoBlocks.ConstGate: Toffoli
using YaoBase: QubitMismatchError
using LinearAlgebra: I

@testset "construction" begin
    @test X + Y + Z == +(X, Y, Z)
    @test sum([X, Y, Z]) == +(X, Y, Z)
    @test X + (Y + Z) == Add(X, Add(Y, Z))

    @test X * Y * Z == *(X, Y, Z)
    @test X * (Y * Z) == *(X, *(Y, Z))
    @test Val(im) * X == Scale(Val(im), X)
    @test X * Val(im) == Scale(Val(im), X)
    @test X * im == Scale(im, X)
    @test im * X == Scale(im, X)
    @test Val(im) * (Val(im) * X) == Scale(Val(-1+0im), X)
    @test im * (im * X) == Scale(-1+0im, X)
    @test (Val(im) * X) * Val(im) == Scale(Val(-1+0im), X)
    @test (im * X) * im == Scale(-1+0im, X)

    @test mat(X - Y) ≈ mat(X) - mat(Y)
    @test mat(X / 2) ≈ 0.5 * mat(X)
end

@testset "Add" begin
    ad = Add(put(3, 2 => X), 3 * put(3, 1 => rot(Y, 0.3)), repeat(3, Z, (2, 3)))
    @test [push!(copy(ad), put(3, 1 => X))...] |> length == 4
    @test insert!(copy(ad), 2, put(3, 1 => X)) isa Add
    @test append!(copy(ad), [put(3, 1 => X)]) |> length == 4
    @test prepend!(copy(ad), [put(3, 1 => X)]) |> length == 4

    @test [occupied_locs(ad)...] |> sort! == [1, 2, 3]
    @test occupied_locs(+(put(5, 2 => X), put(5, 3 => I2))) == (2,)
    ad2 = copy(ad)
    @test ad2 == ad
    push!(ad2, Toffoli)
    @test ad2 != ad
    @test ad2[1:end-1] == ad
    @test ad2[end] == Toffoli
    @test length(ad2) == 4
    println(ad2)
    @test cache_key(ad) != cache_key(ad2)
    @test hash(ad) != hash(ad2)
    ad3 = chsubblocks(ad, [put(3, 1 => X)])
    append!(prepend!(ad3, [rot(Toffoli, 0.2)]), [control(3, 1, 2 => X)])
    @test ad3 == Add(rot(Toffoli, 0.2), put(3, 1 => X), control(3, 1, 2 => X))
    ad3[3] = ad3[1]
    @test ad3 == Add(rot(Toffoli, 0.2), put(3, 1 => X), rot(Toffoli, 0.2))
    @test ad3' == Add(rot(Toffoli, -0.2), put(3, 1 => X), rot(Toffoli, -0.2))

    ad4 = similar(ad3)
    @test typeof(ad4) == typeof(ad3)
    @test length(ad4) == 0

    reg = rand_state(3)
    @test apply!(copy(reg), ad) ≈
          apply!(copy(reg), ad[1]) + apply!(copy(reg), ad[2]) + apply!(copy(reg), ad[3])
    @test mat(ad) * reg.state ≈
          apply!(copy(reg), ad[1]) + apply!(copy(reg), ad[2]) + apply!(copy(reg), ad[3]) |>
          state

    @test Add(3) isa Add
    @test Add(3, [put(3, 3 => X)]) isa Add
    @test_throws QubitMismatchError Add(3, [put(10, 2 => X)])
    @test_throws QubitMismatchError Add(put(10, 2 => X), put(4, 3 => X))
    @test_throws QubitMismatchError apply!(rand_state(2), Add(put(10, 2 => X)))
end

@testset "block operations" begin
    r = rand_state(1)
    @test copy(r) |> X * Y * Z |> state ≈ mat(X * Y * Z) * state(r)
    @test copy(r) |> X + Y + Z |> state ≈ mat(X + Y + Z) * state(r)
end

@testset "empty add" begin
    c = Add(4)
    @test mat(Float64, c) == zeros(16, 16)
end

@testset "algebra" begin
    @test factor(-X) == -1
    @test -(-X) === X
    @test factor(-(2X)) === -2
    @test factor((2X)*(2X)) == 4
    @test factor((2X)*Y) == 2
    @test factor(X*(2Y)) == 2
    @test 2(-X) isa Scale{Int64, 2, Scale{Val{-1}, 2, XGate}}
    @test mat(X^2) ≈ Matrix(I, 2, 2)
end


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
Module for constructing self-contained power system objects.
"""
module PowerSystems

#################################################################################
# Exports

export System
export Topology
export Bus
export Arc
export AggregationTopology
export Area
export LoadZone
export get_aggregation_topology_accessor

export Component
export Device
export get_max_active_power
export get_max_reactive_power
export Branch
export StaticInjection
export ACBranch
export Line
export MonitoredLine
export DCBranch
export HVDCLine
export VSCDCLine
export Transformer2W
export TapTransformer
export PhaseShiftingTransformer

export ThreePartCost
export TwoPartCost
export VariableCost
export MultiStartCost
export MarketBidCost
export get_slopes
export get_breakpoint_upperbounds

export Generator
export HydroGen
export HydroDispatch
export HydroEnergyReservoir
export HydroPumpedStorage

export RenewableGen
export RenewableFix
export RenewableDispatch

export ThermalGen
export ThermalStandard
export ThermalMultiStart

export ElectricLoad
export StaticLoad
export PowerLoad
export PowerLoadPF
export FixedAdmittance
export ControllableLoad
export InterruptibleLoad

export Storage
export GenericBattery
export BatteryEMS

export DynamicComponent
export DynamicInjection
export DynamicGenerator
export DynamicInverter
export DynamicBranch
export HybridSystem
export RegulationDevice

#AVR Exports
export AVR
export AVRFixed
export AVRSimple
export AVRTypeI
export AVRTypeII
export IEEET1
export ESDC1A
export ESDC2A
export ESAC1A
export ESAC6A
export EXAC1
export EXAC1A
export EXAC2
export EXPIC1
export ESST1A
export ESST4B
export SCRX

#Machine Exports
export Machine
export BaseMachine
export RoundRotorMachine
export SalientPoleMachine
export RoundRotorQuadratic
export SalientPoleQuadratic
export RoundRotorExponential
export SalientPoleExponential
export OneDOneQMachine
export MarconatoMachine
export SimpleMarconatoMachine
export AndersonFouadMachine
export SimpleAFMachine
export FullMachine
export SimpleFullMachine

#PSS Exports
export PSS
export PSSFixed
export PSSSimple
export IEEEST

#Shaft Exports
export Shaft
export SingleMass
export FiveMassShaft

#TG Exports
export TurbineGov
export TGFixed
export TGTypeI
export TGTypeII
export GasTG
export GeneralGovModel
export HydroTurbineGov
export IEEETurbineGov1
export SteamTurbineGov1

# Converter Exports
export Converter
export AverageConverter

# DC Source Exports
export DCSource
export FixedDCSource
export ZeroOrderBESS

# Filter Exports
export Filter
export LCLFilter
export LCFilter

# FrequencyEstimator Exports
export FrequencyEstimator
export KauraPLL

# Outer Control Exports
export OuterControl
export VirtualInertia
export ReactivePowerDroop

# InnerControl Export
export InnerControl
export CurrentControl

export Source

export Service
export Reserve
export ReserveDirection
export ReserveUp
export ReserveDown
export ReserveSymmetric
export StaticReserve
export VariableReserve
export AGC
export ReserveDemandCurve
export Transfer
export StaticReserveGroup
export StaticReserveNonSpinning
export VariableReserveNonSpinning

export PTDF
export Ybus
export LODF
export Adjacency
export GeneratorCostModels
export AngleUnits
export BusTypes
export LoadModels
export PrimeMovers
export ThermalFuels
export StateTypes

export TimeSeriesData
export StaticTimeSeries
export Forecast
export AbstractDeterministic
export Deterministic
export Probabilistic
export SingleTimeSeries
export DeterministicSingleTimeSeries
export Scenarios
export ForecastCache
export StaticTimeSeriesCache
export NormalizationFactor
export NormalizationTypes

export get_dynamic_components

export solve_powerflow!
export solve_powerflow
export validate_connectivity

export parse_file
export add_time_series!
export remove_time_series!
export clear_time_series!
export copy_time_series!
export add_component!
export remove_component!
export remove_components!
export clear_components!
export add_service!
export remove_service!
export clear_services!
export get_services
export has_service
export has_time_series
export get_buses
export get_components_in_aggregation_topology
export get_aggregation_topology_mapping
export get_contributing_devices
export get_contributing_device_mapping
export ServiceContributingDevices
export ServiceContributingDevicesKey
export ServiceContributingDevicesMapping
export get_component
export get_components
export get_components_by_name
export get_available_components
export get_forecast_horizon
export get_forecast_initial_timestamp
export get_forecast_interval
export get_forecast_window_count
export get_time_series
export get_time_series_array
export get_time_series_resolution
export get_time_series_timestamps
export get_time_series_values
export get_time_series_names
export get_next_time_series_array!
export get_next_time
export get_horizon
export get_forecast_initial_times
export get_forecast_total_period
export get_resolution
export get_data
export get_lookup
export iterate_components
export get_time_series_multiple
export get_variable_cost
export get_services_bid
export set_variable_cost!
export set_service_bid!
export iterate_windows
export get_window
export transform_single_time_series!
export sanitize_component!
export validate_component
export validate_component_with_system

#export make_time_series
export get_bus_numbers
export get_name
export set_name!
export get_base_power
export get_frequency
export set_units_base_system!
export to_json
export from_json
export serialize
export deserialize
export clear_ext!
export convert_component!
export set_area!
export set_load_zone!
export TamuSystem
export PowerModelsData
export PowerSystemTableData
export add_dyn_injectors!
export get_machine
export get_shaft
export get_avr
export get_prime_mover
export get_pss
export get_converter
export get_outer_control
export get_inner_control
export get_dc_source
export get_freq_estimator
export get_filter
export get_V_ref
export get_P_ref
export get_saturation_coeffs
export set_droop!
export set_participation_factor!
export set_inertia!
export set_reserve_limit_up!
export set_reserve_limit_dn!
export set_cost!
export get_droop
export get_inertia
export get_reserve_limit_up
export get_reserve_limit_dn
export get_participation_factor
export get_cost
export get_units_base
export get_runchecks
export set_runchecks!
export check
export check_component
export check_components

export configure_logging
export open_file_logger
export MultiLogger
export LogEventTracker
export UnitSystem

#################################################################################
# Imports

import SparseArrays
import LinearAlgebra: LAPACK.getri!
import LinearAlgebra: LAPACK.getrf!
import LinearAlgebra: BLAS.gemm
import LinearAlgebra
import Logging
import Dates
import TimeSeries
import DataFrames
import JSON3
import CSV
import YAML
import UUIDs
import Base.to_index
import NLsolve
import InteractiveUtils

import InfrastructureSystems
import InfrastructureSystems:
    Components,
    TimeSeriesData,
    StaticTimeSeries,
    Forecast,
    AbstractDeterministic,
    Deterministic,
    Probabilistic,
    SingleTimeSeries,
    DeterministicSingleTimeSeries,
    Scenarios,
    ForecastCache,
    StaticTimeSeriesCache,
    InfrastructureSystemsComponent,
    InfrastructureSystemsType,
    InfrastructureSystemsInternal,
    DeviceParameter,
    FlattenIteratorWrapper,
    LazyDictFromIterator,
    DataFormatError,
    InvalidRange,
    InvalidValue,
    get_count,
    get_data,
    get_horizon,
    get_resolution,
    get_window,
    get_name,
    set_name!,
    iterate_windows,
    get_time_series,
    get_time_series_array,
    get_time_series_timestamps,
    get_time_series_values,
    get_time_series_names,
    get_next_time_series_array!,
    get_next_time,
    get_units_info,
    set_units_info!,
    to_json,
    from_json,
    serialize,
    deserialize,
    get_time_series_multiple,
    NormalizationFactor,
    NormalizationTypes,
    UnitSystem,
    SystemUnitsSettings,
    open_file_logger,
    validate_struct,
    MultiLogger,
    LogEventTracker

const IS = InfrastructureSystems

#################################################################################

using DocStringExtensions

@template (FUNCTIONS, METHODS) = """
                                 $(TYPEDSIGNATURES)
                                 $(DOCSTRING)
                                 """

#################################################################################
# Includes

"""
Supertype for all PowerSystems components.
All subtypes must include a InfrastructureSystemsInternal member.
Subtypes should call InfrastructureSystemsInternal() by default, but also must
provide a constructor that allows existing values to be deserialized.
"""
abstract type Component <: IS.InfrastructureSystemsComponent end

""" Supertype for "devices" (bus, line, etc.) """
abstract type Device <: Component end

include("definitions.jl")
include("models/static_models.jl")
include("models/dynamic_models.jl")
include("models/injection.jl")

# Include utilities
include("utils/logging.jl")
include("utils/IO/base_checks.jl")

# PowerSystems models
include("models/topological_elements.jl")
include("models/branches.jl")
include("models/operational_cost.jl")
#include("models/network.jl")

# Static types
include("models/services.jl")
include("models/reserves.jl")
include("models/generation.jl")
include("models/storage.jl")
include("models/loads.jl")
include("models/dynamic_generator_components.jl")
include("models/dynamic_inverter_components.jl")
include("models/OuterControl.jl")

# Include all auto-generated structs.
include("models/generated/includes.jl")
include("models/regulation_device.jl")
include("models/HybridSystem.jl")

#Methods for devices
include("models/components.jl")
include("models/devices.jl")

# Dynamic Composed types
include("models/dynamic_generator.jl")
include("models/dynamic_inverter.jl")
include("models/dynamic_machines.jl")
include("models/RoundRotorExponential.jl")
include("models/RoundRotorQuadratic.jl")
include("models/SalientPoleExponential.jl")
include("models/SalientPoleQuadratic.jl")
include("models/dynamic_branch.jl")

include("models/supplemental_constructors.jl")
include("models/supplemental_accessors.jl")

# Definitions of PowerSystem
include("base.jl")
include("data_format_conversions.jl")

#Data Checks
include("utils/IO/system_checks.jl")
include("utils/IO/branchdata_checks.jl")

# cost function TimeSeries convertion
include("models/cost_function_timeseries.jl")

# network calculations
include("utils/network_calculations/common.jl")
include("utils/network_calculations/ybus_calculations.jl")
include("utils/network_calculations/ptdf_calculations.jl")
include("utils/network_calculations/lodf_calculations.jl")

#PowerFlow
include("utils/power_flow.jl")

#Conversions
include("utils/conversion.jl")

# Include Parsing files
include("parsers/common.jl")
include("parsers/enums.jl")
include("parsers/pm_io.jl")
include("parsers/im_io.jl")
include("parsers/power_system_table_data.jl")
include("parsers/power_models_data.jl")
include("parsers/psse_dynamic_data.jl")
include("parsers/TAMU_data.jl")

# Better printing
include("utils/print.jl")

include("models/serialization.jl")

# Download test data
include("utils/data.jl")
import .UtilsData: TestData

end # module


text:
import Base: isvalid

export with_context, set_context!, Context, CONTEXT
export isvalid

abstract type Property end

struct Context
    props::Vector{Property}
end
Context(props::Property...) = Context(collect(props))
Base.vcat(ctx::Context, xs...) = Context(vcat(_props(ctx), _props.(xs)...))
Base.vcat(ctx::Context, ctxs::Context...) = Context(vcat(_props(ctx), _props.(ctxs)...))
_props(ctx::Context) = ctx.props
_props(p) = p
Base.broadcastable(ctx::Context) = Ref(ctx)
set_context!(context::Context) = (global CONTEXT = context)
function with_context(f, context::Context)
    global CONTEXT
    old = CONTEXT
    CONTEXT = context
    try
        f()
    finally
        CONTEXT = old
    end
end

"""
    isvalid(p::Property, ctx::Context) -> Bool

Returns the value corresponding to whether or not the given property object `p` is valid,
given context `ctx`.
"""
isvalid(q::Property, ctx::Context) = any(p -> implies(p, q, ctx), ctx.props)
isvalid(prop::Property) = isvalid(prop, CONTEXT)
isvalid(ctx::Context) = Base.Fix2(isvalid, ctx)

"""
    implies(p::Property, q::Property, , ctx::Context) -> Bool

Determine whether property `p` implies property `q` given context `ctx`.

!!! note
    Since `p` implies `q`, it is implied that `q` is more basic of a property than `p`.
"""
implies(p::Property, q::Property, ctx::Context)

implies(::Property, ::Property, ::Context) = false
implies(p::P, q::P, ::Context) where {P<:Property} = p == q


text:
#
#      Sn - The manifold of the n-dimensional sphere
#  Point is a Point on the n-dimensional sphere.
#
import LinearAlgebra: norm, dot, nullspace
import Base: exp, log, show, cat
export Sphere, SnPoint, SnTVector,show, getValue
export distance, dot, exp, log, manifoldDimension, norm
export randomMPoint, opposite, parallelTransport, zeroTVector
export validateMPoint, validateTVector
export zeroTVector, typeofMPoint, typeofTVector
#
# Type definitions
#
@doc doc"""
    Sphere <: Manifold

The manifold $\mathcal M = \mathbb S^n$ of unit vectors in $\mathbb R^{n+1}$.
This manifold is a matrix manifold (see [`IsMatrixM`](@ref)) and embedded (see
[`IsEmbeddedM`](@ref)).

# Abbreviation

`Sn`

# Constructor

    Sphere(n)

generate the sphere $\mathbb S^n$

`
Its abbreviation is `Sn`.
"""
struct Sphere <: Manifold
  name::String
  dimension::Int
  abbreviation::String
  Sphere(dimension::Int) = new("$dimension-Sphere",dimension,"Sn($(dimension-1))")
end
@doc doc"""
    SnPoint <: MPoint

A point $x$ on the manifold $\mathcal M = \mathbb S^n$ represented by a unit
vector from $\mathbb R^{n+1}$
"""
struct SnPoint{T<:AbstractFloat} <: MPoint
  value::Vector{T}
  SnPoint{T}(value::Vector{T}) where {T <: AbstractFloat} = new(value)
end
SnPoint(value::Vector{T}) where {T <: AbstractFloat} = SnPoint{T}(value)
getValue(x::SnPoint) = x.value;

@doc doc"""
    SnTVector <: TVector

A tangent vector $\xi \in T_x\mathcal M$ on the manifold
$\mathcal M = \mathbb S^n$. For the representation the tangent space can be
given as $T_x\mathbb S^n = \bigl\{\xi \in \mathbb R^{n+1}
\big| \langle x,\xi\rangle = 0\bigr\}$, where $\langle\cdot,\cdot\rangle$
denotes the Euclidean inner product on $\mathbb R^{n+1}$.
"""
struct SnTVector{T <: AbstractFloat} <: TVector
  value::Vector{T}
  SnTVector{T}(value::Vector{T}) where {T <: AbstractFloat} = new(value)
end
SnTVector(value::Vector{T})  where {T <: AbstractFloat}  = SnTVector{T}(value)
getValue(ξ::SnTVector) = ξ.value;
# Traits
# ---
# (a) Sn is a MatrixManifold
@traitimpl IsMatrixM{Sphere}
@traitimpl IsMatrixP{SnPoint}
@traitimpl IsMatrixTV{SnTVector}
# (b) Sn is Embedded
@traitimpl IsEmbeddedM{Sphere}
@traitimpl IsEmbeddedP{SnPoint}
@traitimpl IsEmbeddedV{SnTVector}

# Functions
# ---
@doc doc"""
    distance(M,x,y)

Compute the Riemannian distance on $\mathcal M=\mathbb S^n$ embedded in
$\mathbb R^{n+1}$, which is given by

$ d_{\mathbb S^n}(x,y) = \operatorname{acos} \bigl(\langle x,y\rangle\bigr), $

where $\langle\cdot,\cdot\rangle$ denotes the Euclidean inner product
on $\mathbb R^{n+1}$.
"""
distance(M::Sphere,x::SnPoint,y::SnPoint) = acos(
	min( max(dot(getValue(x), getValue(y) ), -1.), 1. )
)

@doc doc"""
    dot(M,x,ξ,ν)
Compute the Riemannian inner product for two [`SnTVector`](@ref)s `ξ` and `ν`
from $T_x\mathcal M$ of the [`Sphere`](@ref) `M` given by
$\langle \xi, \nu \rangle_x = \langle \xi,\nu \rangle$, i.e. the inner product
in the embedded space $\mathbb R^{n+1}$.
"""
dot(M::Sphere, x::SnPoint, ξ::SnTVector, ν::SnTVector) = dot( getValue(ξ), getValue(ν) )

@doc doc"""
    exp(M,x,ξ[, t=1.0])
Compute the exponential map on the [`Sphere`](@ref) `M`$=\mathbb S^n$ with
respect to the [`SnPoint`](@ref) `x` and the [`SnTVector`](@ref) `ξ`, which can
be shortened with `t` to `tξ`. The formula reads

$\exp_x\xi = \cos(\lVert\xi\rVert_2)x + \sin(\lVert\xi\rVert_2)\frac{\xi}{\lVert\xi\rVert_2}.$
"""
function exp(M::Sphere,x::SnPoint,ξ::SnTVector,t::Float64=1.0)
  len = norm( getValue(ξ) )
  if len < eps(Float64)
    return x
  else
	xV = cos(t*len) * getValue(x)  +  (sin(t*len)/len) * getValue(ξ)
    return SnPoint( xV ./ norm(xV)  )
  end
end
@doc doc"""
    log(M,x,y)
Compute the logarithmic map on the [`Sphere`](@ref)
$\mathcal M=\mathbb S^n$, i.e. the [`SnTVector`](@ref) whose corresponding
[`geodesic`](@ref) starting from [`SnPoint`](@ref) `x` reaches the
[`SnPoint`](@ref)` y` after time 1 on the [`Sphere`](@ref) `M`.
The formula reads for $x\neq -y$

$\log_x y = d_{\mathbb S^n}(x,y)\frac{y-\langle x,y\rangle x}{\lVert y-\langle x,y\rangle x \rVert_2}.$
"""
function log(M::Sphere,x::SnPoint,y::SnPoint)
  scp = dot( getValue(x), getValue(y) )
  ξvalue = getValue(y) - scp*getValue(x)
  ξvnorm = norm(ξvalue)
  if (ξvnorm > eps(Float64))
    return SnTVector( ξvalue * acos(min(max(scp,-1.),1.) )/ξvnorm );
  else
    return zeroTVector(M,x)
  end
end
@doc doc"""
    manifoldDimension(x)
returns the dimension of the [`Sphere`](@ref) `M`$=\mathbb S^n$, the
[`SnPoint`](@ref) `x`, itself embedded in $\mathbb R^{n+1}$, belongs to.
"""
manifoldDimension(x::SnPoint)::Integer = length( getValue(x) )-1
"""
    manifoldDimension(M)
returns the dimension of the [`Sphere`](@ref) `M`.
"""
manifoldDimension(M::Sphere)::Integer = M.dimension
@doc doc"""
    norm(M,x,ξ)
Computes the norm of the [`SnTVector`](@ref) `ξ` in the tangent space
$T_x\mathcal M$ at [`SnPoint`](@ref) `x` of the [`Sphere`](@ref) `M`.
"""
norm(M::Sphere, x::SnPoint, ξ::SnTVector) = norm( getValue(ξ) )
@doc doc"""
    opposite(M,x)
returns the antipodal point of x, i.e. $ y = -x $.
"""
opposite(M::Sphere, x::SnPoint) = SnPoint( -getValue(x) )
@doc doc"""
    parallelTransport(M,x,y,ξ)
Compute the paralllel transport of the [`SnTVector`](@ref) `ξ` from
the tangent space $T_x\mathcal M$ at [`SnPoint`](@ref) `x` to
$T_y\mathcal M$ at [`SnPoint`](@ref)` y` on the [`Sphere`](@ref) `M` provided
that the corresponding [`geodesic`](@ref) $g(\cdot;x,y)$ is unique.
The formula reads

$P_{x\to y}(\xi) = \xi - \frac{\langle \log_xy,\xi\rangle_x}{d^2_{\mathbb S^n}(x,y)}
\bigl(\log_xy + \log_yx \bigr).$
"""
function parallelTransport(M::Sphere, x::SnPoint, y::SnPoint, ξ::SnTVector)
  ν = log(M,x,y)
  νL = norm(M,x,ν)
	return ξ - ( νL > 0 ? dot(M,x,ν,ξ)*(ν + log(M,y,x))/νL^2 : zeroTVector(M,x) )
end
@doc doc"""
    randomMPoint(M [,:Gaussian, σ=1.0])

return a random point on the Sphere by projecting a normal distirbuted vector
from within the embedding to the sphere.
"""
function randomMPoint(M::Sphere, ::Val{:Gaussian}, σ::Float64=1.0)::SnPoint
	v = σ * randn(manifoldDimension(M)+1);
	return SnPoint(v./norm(v))
end
@doc doc"""
    randomTVector(M,x [,:Gaussian,σ=1.0])

return a random tangent vector in the tangent space of the [`SnPoint`](@ref)
`x` on the [`Sphere`](@ref) `M`.
"""
function randomTVector(M::Sphere, x::SnPoint, ::Val{:Gaussian}, σ::Float64=1.0)
    n = σ * randn( size( getValue(x)) ) # Gaussian in embedding
	  nP = n - dot(n,getValue(x))*getValue(x) #project to TpM (keeps Gaussianness)
	  return SnTVector( nP )
end
tangentONB(M::Sphere, x::SnPoint, y::SnPoint) = tangentONB(M,x,log(M,x,y))
function tangentONB(M::Sphere,x::SnPoint,ξ::SnTVector)
    d = manifoldDimension(M)
    A = zeros(d+1,d+1)
    A[1,:] = transpose(getValue(x))
    A[2,:] = transpose(getValue(ξ))
    V = nullspace(A)
    κ = ones(d)
    if ξ != zeroTVector(M,x)
        # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor
		ξ = ξ/norm(M,x,ξ)
		V = cat(getValue(ξ),V,dims=2)
        κ[1] = 0.0 # no curvature along the geodesic direction, if x!=y
    end
    Ξ = [ SnTVector(V[:,i]) for i in 1:d ]
    return Ξ,κ
end
typeofTVector(::Type{SnPoint{T}}) where T = SnTVector{T}
typeofMPoint(::Type{SnTVector{T}}) where T = SnPoint{T}
"""
    typicalDistance(M)

returns the typical distance on the [`Sphere`](@ref)` Sn`: π.
"""
typicalDistance(M::Sphere) = π

@doc doc"""
    validateMPoint(M,x)

validate, whether the [`SnPoint`](@ref) `x` is on the [`Sphere`](@ref) `M`$=\mathbb S^n$,
i.e. that the vector is of the correct dimension $n$ and its norm is $\lVert x \rVert = 1$.
"""
function validateMPoint(M::Sphere, x::SnPoint)
  if length(getValue(x)) ≠ M.dimension+1
    throw( ErrorException(
      "The Point $x is not on the $(M.name), since the vector dimension ($(length(getValue(x)))) is not $(M.dimension+1)."
    ))
  end
  if abs(norm(getValue(x)) - 1) >= 10^(-15)
    throw( ErrorException(
      "The Point $x is not on the $(M.name) since its norm $(norm(getValue(x))) ≠ 1"
    ))
  end
  return true
end

@doc doc"""
    validateTVector(M,x,ξ)

validate, whether the tangent vector [`SnTVector`](@ref) `ξ` is in the tangent
space of [`SnPoint`](@ref) `x` is on the [`Sphere`](@ref) `M`$=\mathbb S^n$,
i.e. that all three lengths are correct and $x^\mathrm{T}\xi = 0$.
"""
function validateTVector(M::Sphere,x::SnPoint,ξ::SnTVector)
  if (length(getValue(x)) ≠ length(getValue(ξ))) || (length(getValue(x)) ≠ M.dimension+1)
    throw( ErrorException(
      "The three dimensions of the $(M.name), the point x ($(length(getValue(x)))), and the tangent vector ($(length(getValue(ξ)))) don't match."
    ))
  end
  if dot(getValue(x),getValue(ξ)) >= 10^(-15)
    throw( ErrorException(
      "The tangent vector $ξ is not a tangent vector to $x on $(M.name) since the inner product is $( dot(getValue(x),getValue(ξ)) ) and not zero."
    ))
  end
  return true
end

@doc doc"""
    ξ = zeroTVector(M,x)

returns a zero vector in the tangent space $T_x\mathcal M$ of the
[`SnPoint`](@ref) $x\in\mathbb S^n$ on the [`Sphere`](@ref)` Sn`.
"""
zeroTVector(M::Sphere, x::SnPoint) = SnTVector(  zero( getValue(x) )  );
# Display
# ---
show(io::IO, M::Sphere) = print(io,"The $(M.name)")
show(io::IO, p::SnPoint) = print(io, "Sn($( getValue(p) ))")
show(io::IO, ξ::SnTVector) = print(io, "SnT($( getValue(ξ) ))")


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    bspdegelev(d::I,c::Matrix{F},k::Vector{F},t::I
              ) where {I<:Integer,F<:AbstractFloat}

Degree elevate a univariate B-Spline.

# Arguments:
- `d`: Degree of the B-Spline.
- `c`: Control points, matrix of size `(dim,nc)`.
- `k`: Knot sequence, vector of size `nk`.
- `t`: Raise the B-Spline degree t times.

# Output:
- `ic`: Control points of the new B-Spline.
- `ik`: Knot vector of the new B-Spline.

# Examples:
```julia
julia> ic,ik = bspdegelev(d,c,k,t)
```

# Reference:
- Algorithm A5.9 from 'The NURBS BOOK' pg206 ([1997Piegl](@cite))
"""
function bspdegelev(d::I,c::Matrix{F},k::Vector{F},t::I
                   ) where {I<:Integer,F<:AbstractFloat}
# TODO: the original C function passed c, k, ic, ik and the point number
#       as reference
# TODO: this function has optimization potential

# added this part as contrary to matlab array sizes are not adaptive
# count unique knots to calculate final array sizes
n0 = k[1];
count = 1;
iknc = 1;
for i in 2:length(k)
  if n0 == k[i]
    iknc += 1;
    continue;
  end # if
  iknc += t+1;
  n0 = k[i];
end # for i
if n0 == k[end];
  iknc += t;
end # if

mc,nc = size(c);

# allocate space for resulting arrays
ic = zeros(F,(mc,iknc-(d+t+1)));
ik = zeros(F,iknc);

n = nc - 1;

bezalfs =  zeros(F,(d+1,d+t+1));
bpts = zeros(F,(mc,d+1));
ebpts = zeros(F,(mc,d+t+1));
Nextbpts = zeros(F,(mc,d+1));
alfs = zeros(F,d);

m = n + d + 1;
ph = d + t;
ph2 = ph ÷ 2;

# compute bezier degree elevation coefficeients
bezalfs[1,1] = 1;
bezalfs[d+1,ph+1] = 1;

for i in 1:ph2
  inv = 1 / bincoeff(ph,i);
  mpi = min(d,i);

  for j in max(0,i-t):mpi
    bezalfs[j+1,i+1] = inv * bincoeff(d,j) * bincoeff(t,i-j);
  end # for j
end # for i

for i in ph2+1:ph-1
  mpi = min(d,i);
  for j in max(0,i-t):mpi
    bezalfs[j+1,i+1] = bezalfs[d-j+1,ph-i+1];
  end # for j
end # for i

mh = ph;
kind = ph + 1;
r = -1;
a = d;
b = d + 1;
cind = 1;
ua = k[1];

for ii in 0:mc-1
  ic[ii+1,1] = c[ii+1,1];
end # for ii

for i in 0:ph
  ik[i+1] = ua;
end # for i

# initialise first bezier seg
for i in 0:d
  for ii in 0:mc-1
    bpts[ii+1,i+1] = c[ii+1,i+1];
  end # for ii
end # for i

# big loop thru knot vector
while b < m
  i = b;
  while b < m && k[b+1] == k[b+2]
    b += 1;
  end # while
  mul = b - i + 1;
  mh +=  mul + t;
  ub = k[b+1];
  oldr = r;
  r = d - mul;

  # insert knot u(b) r times
  if oldr > 0
    lbz = (oldr+2) ÷ 2;
  else
    lbz = 1;
  end # if

  if r > 0
    rbz = ph - (r+1) ÷ 2;
  else
    rbz = ph;
  end # if

  # insert knot to get bezier segment
  if r > 0
    numer = ub - ua;
    for q in d:-1:mul+1
      alfs[q-mul] = numer / (k[a+q+1] - ua);
    end # for q

    for j in 1:r
      save = r - j;
      s = mul + j;

      for q in d:-1:s
        for ii in 0:mc-1
          bpts[ii+1,q+1] = alfs[q-s+1] * bpts[ii+1,q+1] +
                           (1 - alfs[q-s+1]) * bpts[ii+1,q];
        end # for ii
      end # for q

      for ii in 0:mc-1
        Nextbpts[ii+1,save+1] = bpts[ii+1,d+1];
      end # for ii
    end # for j
  end # if
  # end of insert knot

  # degree elevate bezier
  for i in lbz:ph
    for ii in 0:mc-1
      ebpts[ii+1,i+1] = 0;
    end # for ii
    mpi = min(d, i);
    for j in max(0,i-t):mpi
      for ii in 0:mc-1
        ebpts[ii+1,i+1] = ebpts[ii+1,i+1] + bezalfs[j+1,i+1]*bpts[ii+1,j+1];
      end # for ii
    end # for j
  end # for i
  # end of degree elevating bezier

  # must remove knot u=k[a] oldr times
  if oldr > 1
    first = kind - 2;
    last = kind;
    den = ub - ua;
    bet = (ub - ik[kind]) ÷ den;

    # knot removal loop
    for tr in 1:oldr-1
      i = first;
      j = last;
      kj = j - kind + 1;
      # loop and compute the new control points
      # for one removal step
      while j-i > tr
        if i < cind
          alf = (ub - ik[i+1]) / (ua - ik[i+1])
          for ii in 0:mc-1
            ic[ii+1,i+1] = alf * ic[ii+1,i+1] + (1 - alf) * ic[ii+1,i];
          end # for ii
        end # if
        if j >= lbz
          if j-tr <= kind-ph+oldr
            gam = (ub - ik[j-tr+1]) / den;
            for ii in 0:mc-1
              ebpts[ii+1,kj+1] = gam * ebpts[ii+1,kj+1] +
                                 (1 - gam) * ebpts[ii+1,kj+2];
            end # for ii
          else
            for ii in 0:mc-1
              ebpts[ii+1,kj+1] = bet * ebpts[ii+1,kj+1] +
                                 (1 - bet) * ebpts[ii+1,kj+2];
            end # for ii
          end # if
        end # if
        i += 1;
        j -= 1;
        kj -= 1;
      end # while

      first -= 1;
      last += 1;
    end # for tr
  end # if
  # end of removing knot n=k[a]

  # load the knot ua
  if a != d
    for i in 0:ph-oldr-1
      ik[kind+1] = ua;
      kind +=  1;
    end # for i
  end # if

  # load ctrl pts into ic
  for j in lbz:rbz
    for ii in 0:mc-1
      ic[ii+1,cind+1] = ebpts[ii+1,j+1];
    end # for ii
    cind += 1;
  end # for j

  if b < m
  # setup for next pass thru loop
    for j in 0:r-1
      for ii in 0:mc-1
        bpts[ii+1,j+1] = Nextbpts[ii+1,j+1];
      end # for ii
    end # for j
    for j in r:d
      for ii in 0:mc-1
        bpts[ii+1,j+1] = c[ii+1,b-d+j+1];
      end # for ii
    end # for j
    a = b;
    b += 1;
    ua = ub;

  else
  # end knot
    for i in 0:ph
      ik[kind+i+1] = ub;
    end # for i
  end # if
end # while b < m
# End big while loop
return ic,ik;
end # bspdegelev

"""
    bincoeff(n::I,k::I) where{I<:Integer}

Computes the binomial coefficient.

    ( n )      n!
    (   ) = --------
    ( k )   k!(n-k)!

# Arguments:
- `n`: `n` of `n` choose `k`
- `k`: `k` of `n` choose `k`

# Output:
- `b`: binomial coefficient of `n` choose `k`

# Example:
```julia
julia> b = bincoeff(n,k)
```

# Reference:
- Algorithm 6.1.6 from 'Numerical Recipes in C', 2nd Edition pg215
  ([1992Press](@cite))
"""
function bincoeff(n::I,k::I) where{I<:Integer}
# TODO: revisit this later to ommit/replace this with something better?
return floor(0.5 + exp( factln(n) - factln(k) - factln(n-k)));
end # bincoeff

"""
    factln(n::I) where{I<:Integer}

Computes `ln(n!)`.

# Example:
```julia
julia> a = factln(n)
```
"""
function factln(n::I) where{I<:Integer}
# TODO: revisit this later to ommit/replace this with something better?
if n <= 1; return 0; end # if
return log(factorial(n));
end # factln


text:
using MinionsREPL

MinionsREPL.press(831861)


text:
# Superflous code
# Contains code ported from Prof. Hicksn's Mapping_Mod.f90

@doc """
### calcJacobian

!!!!!   NOT NEEDED   !!!!
Calculates the metric Jacobian value based on the B-spline explicit mapping.

**Arguments**

*  `map`  : Object of Mapping type
*  `Jac`  : Array of mesh Jacobians in j,k,m format

"""->

function calcJacobian(map, Jac)

  # Check if map.jkmmax = size(Jac)
  @assert map.jkmmax[1] == size(Jac,1)
  @assert map.jkmmax[2] == size(Jac,2)
  @assert map.jkmmax[3] == size(Jac,3)

  # Allocate and initialize arrays
  scal = zeros(AbstractFloat, 3)
  xi = zeros(AbstractFloat, 3)
  dx = zeros(AbstractFloat, 3, 3)

  scal[:] = 1 ./ (map.jkmmax[:] - 1)

  for m = 1:map.jkmmax[3]
    for k = 1:map.jkmmax[2]
      for j = 1:map.jkmmax[1]
        xi[:] = map.jkmmax[j,k,m,:]
        calcdXdxi(map, xi, jderiv, dx[:,1])
        calcdXdxi(map, xi, jderiv, dx[:,2])
        calcdXdxi(map, xi, jderiv, dx[:,3])
        dx[:,1] *= scal[1]
        dx[:,2] *= scal[2]
        dx[:,3] *= scal[3]
        Jac[j,k,m] = dx[1,1]*dx[2,2]*dx[3,3] + dx[1,2]*dx[2,3]*dx[3,1] +
                     dx[1,3]*dx[2,1]*dx[3,2] - dx[1,1]*dx[2,3]*dx[3,2] -
                     dx[1,2]*dx[2,1]*dx[3,3] - dx[1,3]*dx[2,2]*dx[3,1]
        Jac[j,k,m] = 1/Jac[j,k,m]
      end  # End for j = 1:map.jkmmax[1]
    end  # End for k = 1:map.jkmmax[2]
  end  # End for m = 1:map.jkmmax[3]

  return nothing
end

@doc """
### resizeMapping

!!!! NOT NEEDED FOR FFD  !!!!

Uses edge spacing parameters to refine/coarsen in each parameter direction.
This might be useful for refinement studies where we want to insert or delete
points in a continous way, or for grid sequencing/multigrid applications.

**Arguments**

*  `map`  : Object of Mapping type
*  `jkmmax` : dersired number of (xi, eta, zeta) parameters

"""->

function resizeMapping(map, jkmmax)

  for i = 1:3
    @assert jkmmax[i] > 0
  end

  # Check if resizing is necessary
  if map.jkmmax[1] != jkmmax[1] || map.jkmmax[2] != jkmmax[2] || map.jkmmax[3]
     != jkmmax[3]
    map.xi = zeros(AbstractFloat, jkmmax[1], jkmmax[2], jkmmax[3], 3)
  end

  uvw = zeros(AbstractFloat, 3)

  # Loop over all nodes, calculate each new (xi, eta, zeta)
  for m = 1:map.jkmmax[3]
    uvw[3] = (m-1)/(map.jkmmax[3] - 1)
    for k = 1:map.jkmmax[2]
      uvw[2] = (k-1)/(map.jkmmax[2]-1)
      for j = 1:map.jkmmax[1]

        uvw[1] = (j-1)/(map.jkmmax[1] - 1)
        calcXi(map, uvw, map.xi[j,k,m,:])
        # need to set lower/upper value explicitly to avoid round-off errors
        if j == 1
          map.xi[j,k,m,1] = 0.0
        end
        if k == 1
          map.xi[j,k,m,2] = 0.0
        end
        if m == 1
          map.xi[j,k,m,3] = 0.0
        end
        if j == map.jkmmax[1]
          map.xi[j,k,m,1] = 1.0
        end
        if k == map.jkmmax[2]
          map.xi[j,k,m,2] = 1.0
        end
        if m == map.jkmmax[3]
          map.xi[j,k,m,3] = 1.0
        end

      end  # End for j = 1:map.jkmmax[1]
    end  # End for k = 1:map.jkmmax[2]
  end  # End for m = 1:map.jkmmax[3]

  return nothing
end  # End function resizeMapping(map, jkmmax)

@doc """
### calcXi

!!!! NOT NEEDED FOR FFD !!!

Calculates the intermediate parameters, (xi,eta,zeta), using the edge spacing
control functions defined by map%edgpar(1:2,:,:)

**Arguments**

*  `map`  : Object of Mapping type
*  `uvw`  : j,k,m indices scaled to fit in [0,1],  i.e. uniform spacing
*  `xi`   : Parametric coordinates (length = 3)

"""->

function calcXi(map, uvw, xi)

  # Allocate and initialize arrays
  # xi = zeros(AbstractFloat, 3)
  xiedg = zeros(AbstractFloat, 4, 3)
  F = zeros(AbstractFloat, 3)
  dF = zeros(AbstractFloat, 3, 3)
  dxi = zeros(AbstractFloat, 3)
  psi = zeros(AbstractFloat, 3)
  dpsi = zeros(AbstractFloat, 3)

  for di = 1:3
    it1 = mod(di,3) + 1
    it2 = mod(di+1,3) + 1
    for edg 1:4
      A = map.edge_param[1, edg, di]
      b = map.edge_param[2, edg, di]
      fac = uvw[di] - 0.5
      if b == 0.0
        # Use taylor series to determine u
        fac = 1.0 + 2.0*fac
      else
        fac = 1.0 + tanh(b*fac)/tanh(0.5*b)
      end  # End if b == 0.0
      xiedg[edg, di] = fac/( 2*A + (1 - A)*fac )
    end  # End for edg 1:4
  end  # End for di = 1:3

  # calculate the intermediate params by find the intersection of 3 trilinear
  # surface (via Newton's method)

  # Get initial estimate
  for di = 1:3
    xi[di] = 0.0
    for edg = 1:4
      xi[di] = 0.25*xiedg[edg,di]
    end
  end

  # loop until sufficiently converged

  for n = 1:100
    # Get the RHS and LHS for Newton's method
    F[:] = -xi[:]
    fill!(dF, 0.0)
    for di = 1:3
      it1 = mod(di,3) + 1
      it2 = mod(di+1,3) + 1
      dF[di,di] = 1.0
      for edg = 1:4
        if mod(edg,2) == 1
          psi[it1] = 1.0 - xi[it1]
          dpsi[it1] = 1.0
        else
          psi[it1] = xi[it1]
          dpsi[it1] = 1.0
        end  # End if mod(edg,2) == 1

        if mod((edg+1)/2,2) == 1
          psi[it2] = 1.0 - xi[it2]
          dpsi[it2] = -1.0
        else
          psi[it2] = xi[it2]
          dpsi[it2] = 1.0
        end  # End if mod((edg+1)/2,2) == 1

        F[di] += psi[it1]*psi[it2]*xiedg[edg,di]
        dF[di,it1] -= -dpsi[it1]*psi[it2]*xiedg[edg,di]
        dF[di,it2] -= -psi[it1]*dpsi[it2]*xiedg[edg,di]
      end  # End for edg = 1:4
    end # End for di = 1:3

    # Solve for the update
    # TODO: Solve for the update

    # update parameter and exit loop
    xi[:] += F[:]
    if norm(F,2) < 1e-13
      return nothing
    end

  end  # End for n = 1:100

  return nothing
end

# parameter functions
@doc """
### calcChordLengthParam

!!! NOT NEEDED FOR FFD !!!

Calculates the mapping parameter values based on the chord lengths in physical
xyz space. It also calculates the edge knot vectors.

**Inputs**

*  `map` : Object of Mapping type
*  `xyz` : Mesh in physical space on which the chord length parameterization is
           based on

**Outputs**

*  None

""" ->

function calcChordLengthParam(map, xyz)

  # Intermediate variables
  jkm = zeros(Int, 3)
  vec = zeros(3)

  # Loop over each direction, Find the parameter values for a curve of constant
  # it1 & it2
  for di = 1:3
    it1 = mod(di,3) + 1
    it2 = mod(di+1,3) + 1

    for jit1 = 1:jkmmax[it1]
      jkm[it1] = jit1
      for jit2 = 1:jkmmax[it2]
        jkm[it2] = jit2
        jkm[di] = 1
        map.xi[jkm[1], jkm[2], jkm[3], di] = 0.0
        for jdi = 1:(map.jkmmax[di]-1)
          jkm[di] = jdi
          j = jkm[1]
          k = jkm[2]
          m = jkm[3]
          jkm[di] = jdi + 1
          jp = jkm[1]
          kp = jkm[2]
          mp = jkm[3]

          # Get the arc-length element
          vec[:] = xyz[jp, kp, mp, :] - xyz[j, k, m, :]
          ds = norm(vec,2)
          map.xi[jp,kp,mp,di] = map.xi[j,k,m,di] + ds
        end  # End for jdi = 1:(map.jkmmax[di]-1)

        # Normalize the arc-length parameters
        jkm[di] = map.jkmmax[di]
        j = jkm[1]
        k = jkm[2]
        m = jkm[3]
        fac = 1 / map.xi[j,k,m,di]
        map.xi[j,k,m,di] = 1.0
        for jdi = 2:(map.jkmmax[di]-1)
          jkm[di] = jdi
          j = jkm[1]
          k = jkm[2]
          m = jkm[3]
          map.xi[j,k,m,di] = fac*map.xi[j,k,m,di]
        end  # End for jdi = 2:(map.jkmmax[di]-1)
      end  # End for jit2 = 1:jkmmax[it2]
    end  # End for jit1 = 1:jkmmax[it1]
  end  # end for di = 1:3

  # At this point, we have arc length parameters, so arc-length-based edge knots
  # can be found

  for di = 1:3
    it1 = mod(di,3) + 1
    it2 = mod(di+1,3) + 1
    order = map.order[di]
    nctl = map.nctl[di]

    for edg = 1:4
      if mod(edg,2) == 1
        jkm[it1] = 1
      else
        jkm[it1] = map.jkmmax[it1]
      end  # End if mod(edg,2) == 1

      if mod((edg+1)/2,2) == 1
        jkm[it2] = 1
      else
        jkm[it2] = map.jkmmax[it2]
      end  # End if mod((edg+1)/2,2) == 1

      # Set the edge knot vector's end knots
      map.edge_knot[1:order, edg, di] = 0.0
      map.edge_knot[nctl+1:nctl+order, edg, di] = 1.0

      # Find point intervals per knot intervals
      ppk = (map.jkmmax[di]-1) / (nctl+1-order)

      for i = order+1:nctl
        knotloc = ppk*(i-order)  # knot location
        for jdi = 1:(map.jkmmax[di]-1)
          jkm[di] = jdi

          if knotloc >= jdi-1 && knotloc < jdi
            jkm[di] = jdi
            j = jkm[1]
            k = jkm[2]
            m = jkm[3]
            jkm[di] = jdi + 1
            jp = jkm[1]
            kp = jkm[2]
            mp = jkm[3]
            map.edge_knot[i,edg,di] = (map.xi[jp,kp,mp,di] - map.xi[j,k,m,di])*
                                      (knotloc - jdi + 1) + map.xi[j,k,m,di]
          end  # End if knotloc >= jdi-1 && knotloc < jdi
        end # End for jdi = = 1:(map.jkmmax[di]-1)
      end  # End for i = order+1:nctl
    end  # End for edg = 1:4
  end # End for di = 1:3

  return nothing
end  # end function calcChordLengthParam

@doc """
### calcEdgeSpacingParam

!!!  NOT NEEDED  !!!
Uses the xi, eta, zeta values to calculate the parameters for hyperbolic tangent
spacing function. These parameters are

A = sqrt(sp2/sp1)
b, where, sinh(b) - b/((N-1) \* sp1 \* sp2) = 0

and sp1 and sp2 are the (xi, eta, zeta) spacing at the ends of the edges.
A is stored in `map.edge_param[1,:,:]`, and b in `map.edge_param[2, :, :]`

**Arguments**

*  `map` : Object of Mapping type

""" ->

function calcEdgeSpacingParam(map)

  # Initialize working array
  jkm = zeros(Int,3)

  for di = 1:3
    it1 = mod(di,3) + 1
    it2 = mod(di+1,3) + 1

    for edg = 1:4
      if mod(edg,2) == 1
        jkm[it1] = 1
      else
        jkm[it1] = map.jkmmax[it1]
      end

      if mod((edg+1)/2,2) == 1
        jkm[it2] = 1
      else
        jkm[it2] = map.jkmmax[it2]
      end

      # Find the spacing at the ends of the edge edg
      jkm[di] = 1
      j = jkm[1]
      k = jkm[2]
      m = jkm[3]
      jkm[di] = 2
      jp = jkm[1]
      kp = jkm[2]
      mp = jkm[3]
      dx2 = map.xi[jp,kp,mp,di] - map.xi[j,k,m,di]

      ### #if 0 satement is supposed to come here
      ###
      ###

      # Set parameter A
      map.edge_param[1, edg,di] = sqrt(dx2/dx1)

      # Set parameter b
      Nm1 = map.jkmmax[di] - 1
      x2 = (map.jkmmax[di] - 1)*sqrt(dx1*dx2)
      if x2 > 1
        # 1/(N-1) is smaller than sqrt(dx1*dx2), so b will be zero
        map.edge_param[2, edg, di] = 0.0
      else
        # need to find b using Newton's method
        x1 = small
        x2 = sqrt(6*(1/x2 - 1))
        map.edge_param[2,edg,di] =  bfuncNewtonSolve(x2, Nm1, dx1, dx2)
      end  # End if x2 > 1
    end    # End for edg = 1:4
  end      # End for di = 1:3

  return nothing
end  # End function calcEdgeSpacingParam(map)

function bfunc(x, Nm1, dx1, dx2)

  tmp = 1 / (Nm1*sqrt(dx1*dx2))
  f = sinh(x) - x*tmp
  df = cosh(x) - tmp

  return f, df
end  # end function bfunc

function bfuncNewtonSolve(guess, Nm1, dx1, dx2)

  b = guess
  b_next = zero(AbstractFloat)

  for i = 1:50
    f, df = bfunc(b, Nm1, dx1, dx2)
    b_next = b - f/df
    if norm(b_next-b, 2) < 1e-12
      break
    else
      b = b_next
    end
  end # End for i = 1:50

  return b_next
end

# Grid.jl

####  -------------------

#  NOT NEEDED !!!!!!!!!!

#### ---------------------

@doc """
### calcGrid

Calculates the mesh coordinates based on the element node values in the map,
hermite basis functions and the refinement level. The basis function values of
map should be known

**Arguments**

*  `map` : Object of Mapping type
*  `xyz` : Array of mesh node coordinates in j,k,m,1:3 format

"""->

function calcGrid(map, xyz)

  @assert size(map.jkmmax) == size(xyz)

  for m = 1:map.jkmmax[3]
    for k = 1:map.jkmmax[2]
      for j = 1:map.jkmmax[1]
        xi =sview(map.ki, j, k, m, :)
        dXdxi(map, xi, jderiv, xyz[j,k,m,:])
      end  # End for j = 1:map.jkmmax[1]
    end  # End for k = 1:map.jkmmax[2]
  end  # End for m = 1:map.jkmmax[3]

  return nothing
end

@doc """
### contractWithdGdB

Used to contract the matrix dG/dB, the derivative of the mesh node coordinates
with respect to the mapping control points, with a given array, here labelled
dJdGrid and stored in (j,k,m,:) format. This can be used to find
dJ/dB = (dJ/dG)(dG/dB) for some objective J. In practice, the objective used is
actually

             J_practice = J_obj + psi^T dot R

where J_obj is the actual objective, psi are the flow adjoints, and R is the
flow residual.

**Arguments**

*  `map`     : Object of mapping type
*  `dJdGrid` : The derivative of the objective w.r.t the grid coordintes

"""->

function contractWithdGdB(map, dJdGrid)

  @assert map.numnodes[1] == size(dJdGrid, 1)
  @assert map.numnodes[2] == size(dJdGrid, 2)
  @assert map.numnodes[3] == size(dJdGrid, 3)

  # Allocate and initialize arrays
  basis = zeros(AbstractFloat, maximum(map.order), 3)
  fill!(map.work, 0.0)

  # Loop over the nodes of the mapping
  for m = 1:map.numnodes[3]
    for k = 1:map.numnodes[2]
      for j = 1:map.numnodes[1]

        # Store the parameter values and dJdGrid
        xi[:] = map.xi[j,k,m,:]
        dJdG =sview(dJdGrid, j, k, m, 1:3)

        # Evaluate the knot vectors and basis values
        # TODO: Compute the knot vector for subsequent operations
        for di = 1:3
          span[di] = findSpan(xi[di], map, di)
          basisFunctions(map, basis[:,di], di, xi[di], span[di])
        end  # End for di = 1:3

        for r = 1:map.order[3]
          rs = r + span[3] - map.order[3]
          for q = 1:map.order[2]
            qs = q + span[2] - map.order[2]
            for p = 1:map.order[1]

              ps = p + span[1] - map.order[1]
              coeff = basis[p,1]*basis[q,2]*basis[r,3]
              map.work[ps, qs, rs, 1:3] += coeff*dJdG[:]

            end # End for p = 1:map.order[1]
          end   # End for q = 1:map.order[2]
        end     # End for r = 1:map.order[3]

      end  # End for j = 1:map.numnodes[1]
    end    # End for k = 1:map.numnodes[2]
  end      # End for m = 1:map.numnodes[3]

  return nothing
end  # End function contractWithdGdB(map, dJdGrid)

@doc """
### fitGrid

Finds the control point values such that, when evaluated at in computational
space, the mapping produces the mesh xyz as close as possible. Parameter
correction was added by L. Olague and J. Hicken and ported to the language
`Julia`

**Arguments**

*  `map` : Object of Mapping type
*  `xyz` : Mesh to be approximated. On exit, htis is over written with the new
           mesh produced by the map
*  `L2_error` : The L2 error squared in fitting the grid is added to this
                variable at the end of routine
*  `max_error`: The values in this array are updated if the infinity norm in the
                fit of this grid is larger

"""

function fitGrid(map, xyz, L2_error, max_error)

  # Allocate and initialize arrays

  return nothing
end


text:
if Winston.output_surface == :tk
    using Tk
elseif Winston.output_surface == :gtk
    using Gtk
    include("canvasgtkcompat.jl")
else
    error("unsupported output_surface for this test")
end

win = Toplevel("TestCanvas", 400, 200)
c = Canvas(win, 400, 200)
pack(c, expand=true, fill="both")

p = FramedPlot()
x = range(0,stop=10,length=1001)
y = sin.(x)
add(p, Curve(x, y, "color", colorant"green"))

Winston.display(c, p)


text:
################################################################################
# Active Set
################################################################################

function active(game_con::GameConstraintValues, stamp::CStamp)
	probsize = game_con.probsize
    N = probsize.N
    p = probsize.p
    px = probsize.px

    a = BitArray([0])
    if stamp.con == :col
        if valid(stamp, N, p)
            i = stamp.i
            j = stamp.j
            j_ = findfirst(x -> all(x.con.x2 .== px[j]), game_con.state_conval[i])
            conval = game_con.state_conval[i][j_]
            l = findfirst(x -> x==stamp.k, conval.inds)
            a[1] = conval.active[l][1]
			# if conval.vals[l][1] >= 0. || conval.λ[l][1] != 0.
                # a[1] = 1
            # end
        end
    end
    return a
end

function active_vertical_mask!(ascore::ActiveSetCore, game_con::GameConstraintValues)
    probsize = game_con.probsize
    N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	S = probsize.S

    ascore.vmask = Vector{Int}(1:S)
	stamp = CStamp()
	for k = 2:N
		for i = 1:p
			for j = i+1:p
				stampify!(stamp, :v, :col, i, j, k)
				if active(game_con, stamp)[1] == 1
					ascore.vmask = vcat(ascore.vmask, vertical_idx(ascore, stamp))
				end
			end
		end
	end
    return nothing
end

function active_horizontal_mask!(ascore::ActiveSetCore, game_con::GameConstraintValues)
    probsize = game_con.probsize
    N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	S = probsize.S

    ascore.hmask = Vector{Int}(1:S)
	stamp = CStamp()
	for k = 2:N
		for i = 1:p
			for j ∈ setdiff(1:p,i)
				stampify!(stamp, :h, :col, i, j, k)
				if active(game_con, stamp)[1] == 1
					ascore.hmask = vcat(ascore.hmask, horizontal_idx(ascore, stamp))
				end
			end
		end
	end
    return nothing
end

################################################################################
# Residual
################################################################################

function get_collision_conval(game_con::GameConstraintValues, i::Int, j::Int)
	px = game_con.probsize.px
	j_ = findfirst(
		x -> typeof(x.con) <: CollisionConstraint && all(x.con.x2 .== px[j]),
		game_con.state_conval[i]
	)
	if j_ == nothing
		valid = false
		return nothing, valid
	else
		valid = true
		return game_con.state_conval[i][j_], valid
	end
end

function residual!(ascore::ActiveSetCore, prob::GameProblem{KN,n,m,T,SVd,SVx},
	pdtraj::PrimalDualTraj{KN,n,m,T,SVd}) where {KN,n,m,T,SVd,SVx}

	probsize = ascore.probsize
	N = probsize.N
	p = probsize.p
	px = probsize.px
	S = probsize.S

	# Initialization
	ascore.res .= 0.0

	residual!(prob, pdtraj)
	ascore.res[1:S] .= prob.core.res

	evaluate!(prob.game_con, pdtraj.pr)
	stamp = CStamp()
	for i = 1:p
		for j = i+1:p
			conval, v = get_collision_conval(prob.game_con, i, j)
			if v
				for (l,k) in enumerate(conval.inds)
					stampify!(stamp, :v, :col, i, j, k)
					valid(stamp,N,p) ? add2sub(ascore.res_sub[stamp], conval.vals[l]) : nothing
				end
			end
		end
	end
    return nothing
end

################################################################################
# Residual Jacobian
################################################################################

function residual_jacobian!(ascore::ActiveSetCore, prob::GameProblem{KN,n,m,T,SVd,SVx},
	pdtraj::PrimalDualTraj{KN,n,m,T,SVd}) where {KN,n,m,T,SVd,SVx}

	probsize = ascore.probsize
	N = probsize.N
	p = probsize.p
	S = probsize.S

	# Reset!
	sparse_zero!(ascore.jac)

	residual_jacobian!(prob, pdtraj)
	ascore.jac[1:S,1:S] .= prob.core.jac

    # Allocations
	vs = VStamp()
	hs = HStamp()
    cs = CStamp()

	jacobian!(prob.game_con, pdtraj.pr)
	for i = 1:p
		for j ∈ setdiff(1:p,i)
			conval, v = get_collision_conval(prob.game_con, i, j)
			if v
				for (l,k) in enumerate(conval.inds)
					stampify!(vs, :opt, i, :x, 1, k)
					stampify!(cs, :h, :col, i, j, k)
					if valid(cs,N,p) && valid(vs,N,p)
						add2sub(ascore.jac_sub[(vs,cs)], conval.jac[l]')
					end
					if valid(hs,N,p) && valid(cs,N,p)
						stampify!(hs, :x, 1, k)
						stampify!(cs, :v, :col, i, j, k)
						add2sub(ascore.jac_sub[(hs,cs)], conval.jac[l])
					end
				end
			end
		end
	end
    return nothing
end

function update_nullspace!(ascore::ActiveSetCore, prob::GameProblem{KN,n,m,T,SVd,SVx},
	pdtraj::PrimalDualTraj{KN,n,m,T,SVd}; atol::T=1e-20) where {KN,n,m,T,SVd,SVx}

	update_active_set!(prob.game_con, pdtraj.pr)
	active_vertical_mask!(ascore, prob.game_con)
	active_horizontal_mask!(ascore, prob.game_con)
	residual_jacobian!(ascore, prob, pdtraj)
	djac = Matrix(ascore.jac[ascore.vmask, ascore.hmask])
	reset!(ascore.null)
	add_matrix!(ascore.null, nullspace(djac, atol=1e-20), ascore.hmask)
	return nothing
end


text:
module TCX
using EzXML, Dates, DataFrames, Geodesy, Mocking
import Base.show

export parse_tcx_dir, parse_tcx_file, getActivityType, getDataFrame, getDistance, getDistance2, getDuration, getAverageSpeed, getAveragePace

const OK = 200
const CLIENT_ERROR = 400
const CLIENT_TCX_ERROR = 401
const NOT_FOUND = 404
const SERVER_ERROR = 500

struct TrackPoint
    Time::DateTime
    Latitude::Float64
    Longtitude::Float64
    HeartRateBpm::Int32
    AltitudeMeter::Float64
    DistanceMeter::Float64
end

struct TCXRecord
    Id::DateTime
    Name::String
    ActivityType::String
    DistanceStatic::Float64
    DurationStatic::Float64
    HeartRate::Int32
    TrackPoints::Array{TrackPoint}
end

function parse_tcx(tcxdoc::EzXML.Document)
    root_element = root(tcxdoc)
    # Check if TCX
    if nodename(root_element) != "TrainingCenterDatabase"
        return CLIENT_TCX_ERROR, nothing
    end

    # Type - "/*/*[1]/*[1]/@Sport"
    aType = nodecontent(findfirst("/*/*[1]/*[1]/@Sport", tcxdoc))
    # Id - "/*/*[1]/*/*[1]"
    xid = nodecontent(findfirst("/*/*[1]/*/*[1]", tcxdoc))

    aId = convertToDateTime(xid)
    # Name = "/*/*[1]/*[1]/*[3]"
    aName = nodecontent(findfirst("/*/*[1]/*/*[2]", tcxdoc))
    # Lap - "/*/*[1]/*/*[2]"
    # TotalSeconds - "/*/*[1]/*/*[2]/*[1]"
    aTime = parse(Float64, nodecontent(findfirst("/*/*[1]/*/*[2]/*[1]", tcxdoc)))
    aDistance = parse(Float64, nodecontent(findfirst("/*/*[1]/*/*[2]/*[2]", tcxdoc)))
    # DistanceMeters - "/*/*[1]/*/*[2]/*[2]"
    # AverageHeartRateBpm - "/*/*[1]/*/*[2]/*[5]/*[1]"
    xbpm = findfirst("/*/*[1]/*/*[2]/*[5]/*[1]", tcxdoc)
    if xbpm === nothing
        aHeartRateBpm = 0
    else
        aHeartRateBpm = parse(Int32, nodecontent(xbpm))
    end
    # TrackPoints - "/*/*[1]/*/*[2]/*[9]/*"
    tp_Points = findall("/*/*[1]/*/*[2]/*[9]/*", tcxdoc)
    aTrackPoints = Array{TrackPoint, size(tp_Points, 1)}[]
    for tp in tp_Points
        xtime = nodecontent(findfirst("./*[local-name()='Time']", tp))
        tp_time = convertToDateTime(xtime)
        tp_lat = parseNode(Float64, "./*[local-name()='Position']/*[local-name()='LatitudeDegrees']", tp)
        tp_lont = parseNode(Float64, "./*[local-name()='Position']/*[local-name()='LatitudeDegrees']", tp)
        tp_bpm = parseNode(Int32, "./*[local-name()='HeartRateBpm']/*[1]", tp)
        tp_dist = parseNode(Float64, "./*[local-name()='TPX']", tp)
        tp_alt = parseNode(Float64, "./*[local-name()='AltitudeMeters']", tp)

        aTrackPoints = vcat(aTrackPoints, TrackPoint(tp_time, tp_lat, tp_lont, tp_bpm, tp_dist, tp_alt))
    end

    return OK, TCXRecord(aId, aName, aType, aDistance, aTime, aHeartRateBpm, aTrackPoints)
end

function parse_tcx_str(str::String)
    try
        status, parsed_tcx = parse_tcx(EzXML.parsexml(str))
        warn_on_tcx_error(status, str, false)
        return status, parsed_tcx
    catch e
        if isa(e, EzXML.XMLError)
            @error "Invalid XML string: $str"
            return CLIENT_ERROR, nothing
        end
    end
end

function parse_tcx_file(file::String)
    file_path = abspath(file)
    if isfile(file_path) == false
        return NOT_FOUND, nothing
    end
    xmldoc = try @mock EzXML.readxml(file_path)
    catch e
       if isa(e, EzXML.XMLError)
           # Not a valid XML document
           @warn "Invalid XML document: $file_path"
           return CLIENT_ERROR, nothing
       else
           throw(e)
       end
    end

    status, parsed_tcx = parse_tcx(xmldoc)
    warn_on_tcx_error(status, file_path, true)

    return status, parsed_tcx
end

function warn_on_tcx_error(status::Int, thing::String, isFile::Bool)
    if status == CLIENT_TCX_ERROR
        @warn "Invalid TCX $(isFile ? document : string): $(thing)"
    end
end

#=
= Parses an XML node based on an XPATH and data type.
=
= If the node is a `nothing` value, the function returns the data
= type's version of 0
=#
function parseNode(dType, path, node)
    node_check = findfirst(path, node)
    if node_check !== nothing
        return parse(dType, nodecontent(node_check))
    else
        return dType(0)
    end
end

function parse_tcx_dir(path::String)
    if ispath(path) == false
        @warn "Invalid path: $path"
        return SERVER_ERROR, nothing
    end

    tcxArray = Array{TCXRecord}[]
    searchdir(path, key) = filter(x->occursin(key, x), readdir(path))

    for f in searchdir(path, ".tcx")
        err, tcx = parse_tcx_file(joinpath(path, f))
        if err == OK
            tcxArray = vcat(tcxArray, tcx)
        end
    end

    if length(tcxArray) > 0
        return OK, tcxArray
    else
        return NOT_FOUND, nothing
    end
end

function getActivityType(record::TCXRecord)
    return record.ActivityType
end

function getDataFrame(record::TCXRecord)
    return DataFrame(record.TrackPoints)
end

function getDataFrame(tcxArray::Array{Any, 1})
    aTP = Array{TrackPoint}[]
    for t in tcxArray
        aTP = vcat(aTP, t.TrackPoints)
    end
    return DataFrame(aTP)
end

function getDistance(record::TCXRecord)
    return record.DistanceStatic
end

function getDistance2(record::TCXRecord)
    total_distance = 0
    df = getDataFrame(record)
    num_of_rows = size(df, 1)
    for i in 1:num_of_rows
        if i < num_of_rows
            total_distance += distance(
                                       LLA(df[i, :Latitude], df[i, :Longtitude], df[i, :AltitudeMeter]),
                                       LLA(df[i+1, :Latitude], df[i+1, :Longtitude], df[i+1, :AltitudeMeter])
               )
        end
    end
    return total_distance
end

function getAverageSpeed(record::TCXRecord)
    return (record.DistanceStatic /1000) / (record.DurationStatic / 3600)  # km/h
end

function getAveragePace(record::TCXRecord)
    return (record.DurationStatic / 60) / (record.DistanceStatic / 1000) # min/km
end

function getDuration(record::TCXRecord)
    return record.DurationStatic
end

#=
= Converts a datetime string into the proper datetime based on string length.
=
= Will assume that an ArgumentError is due to
= https://github.com/JuliaLang/julia/issues/23049 and will attempt to work
= around this.
=#
function convertToDateTime(datestr::String)::DateTime
    m = match(r"\d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}(Z?|\.\d{1,3}Z?)", datestr)
    format_prefix = "yyyy-mm-ddTHH:MM:SS"
    if m === nothing
        msg = "'$(datestr)' is improperly formatted. Must be in the form "
        msg = msg * "'$(format_prefix)Z' or '$(format_prefix).sssZ'"
        throw(ArgumentError(msg))
    else
        suffix = replace(m.captures[1], r"\d" => "s")
        try
            return @mock DateTime(m.match, format_prefix * suffix)
        catch e
            if isa(e, ArgumentError)
                # OK! FINE! NO Z THEN!
                return DateTime(m.match[1:end-1], format_prefix * (suffix[1:end-1]))
            else
                throw(e)
            end
        end
    end
end

Base.show(io::IO, tcx::TCXRecord) = print(io, "$(tcx.ActivityType) $(tcx.DistanceStatic/1000) km at $(tcx.Id) for $(tcx.DurationStatic) seconds.")
end #module_end


text:
using SossModels
using Test

@testset "SossModels.jl" begin
    # Write your own tests here.
end


text:
@testset "Property type" begin
    atom = Atom("")
    set_property!(atom, "test", false)
    @test property(atom, "test") == false

    frame = Frame()
    set_property!(frame, "test", 1234.567)
    @test property(frame, "test") == 1234.567

    set_property!(atom, "test2", "TESTINGTESTING")
    @test property(atom, "test2") == "TESTINGTESTING"

    set_property!(frame, "test2", [1.0,2.0,3.0])
    @test property(frame, "test2") == [1.0,2.0,3.0]

    lipsum = "Enim possimus non tenetur laborum. Quo ut
    reiciendis in at suscipit in nulla quas. Eos minus animi quis quo enim
    facilis dolores quo.

    Iusto incidunt molestiae distinctio cum officia cupiditate tempore qui. Sunt
    error voluptatibus reprehenderit rerum aut a. Earum illum asperiores est
    accusamus illo sit.

    Magni incidunt enim nam. Nisi adipisci harum corrupti reprehenderit nihil.
    Voluptatem aut non labore et sint nesciunt qui. Perspiciatis atque laborum
    aut."
    set_property!(atom, "very long", lipsum)
    @test property(atom, "very long") == lipsum
end


text:
@with_kw mutable struct DRegularizedData{T <: Real}
    Q::T = 1e10
    Q̃::T = 1e10
    θ::T = -1e10
    σ::T = 1.0
    timestamp::Int = 1
    iterations::Int = 1
    major_iterations::Int = 0
    minor_iterations::Int = 0
end

@with_kw mutable struct DRegularizedParameters{T <: Real}
    κ::T = 0.3
    τ::T = 1e-6
    σ::T = 1.0
    σ̅::T = 4.0
    σ̲::T = 0.5
    log::Bool = true
    autotune::Bool = false
end

struct DRegularized{T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver} <: AbstractLShapedSolver{T,A,M,S}
    structuredmodel::JuMP.Model
    solverdata::DRegularizedData{T}

    # Master
    mastersolver::M
    mastervector::A
    c::A
    x::A
    Q_history::A

    committee::Vector{SparseHyperPlane{T}}
    inactive::Vector{SparseHyperPlane{T}}
    violating::PriorityQueue{SparseHyperPlane{T},T}

    # Subproblems
    nscenarios::Int
    subobjectives::Vector{A}
    finished::Vector{Int}

    # Workers
    subworkers::Vector{SubWorker{T,A,S}}
    work::Vector{Work}
    decisions::Decisions{A}
    cutqueue::CutQueue{T}

    # Trust region
    ξ::A
    Q̃_history::A
    σ_history::A

    # Cuts
    θs::A
    cuts::Vector{SparseHyperPlane{T}}
    θ_history::A

    # Params
    parameters::DRegularizedParameters{T}
    progress::ProgressThresh{T}

    @implement_trait DRegularized IsRegularized
    @implement_trait DRegularized IsParallel

    function (::Type{DRegularized})(model::JuMP.Model,ξ₀::AbstractVector,mastersolver::AbstractMathProgSolver,subsolver::AbstractMathProgSolver; kw...)
        if nworkers() == 1
            warn("There are no worker processes, defaulting to serial version of algorithm")
            return Regularized(model,ξ₀,mastersolver,subsolver; kw...)
        end
        length(ξ₀) != model.numCols && error("Incorrect length of starting guess, has ",length(ξ₀)," should be ",model.numCols)
        !haskey(model.ext,:SP) && error("The provided model is not structured")

        T = promote_type(eltype(ξ₀),Float32)
        c_ = convert(AbstractVector{T},JuMP.prepAffObjective(model))
        c_ *= model.objSense == :Min ? 1 : -1
        mastervector = convert(AbstractVector{T},copy(ξ₀))
        x₀_ = convert(AbstractVector{T},copy(ξ₀))
        ξ₀_ = convert(AbstractVector{T},copy(ξ₀))
        A = typeof(x₀_)

        msolver = LQSolver(model,mastersolver)
        M = typeof(msolver)
        S = LQSolver{typeof(LinearQuadraticModel(subsolver)),typeof(subsolver)}
        n = StochasticPrograms.nscenarios(model)

        lshaped = new{T,A,M,S}(model,
                               DRegularizedData{T}(),
                               msolver,
                               mastervector,
                               c_,
                               x₀_,
                               A(),
                               convert(Vector{SparseHyperPlane{T}},linearconstraints(model)),
                               Vector{SparseHyperPlane{T}}(),
                               PriorityQueue{SparseHyperPlane{T},T}(Reverse),
                               n,
                               Vector{A}(),
                               Vector{Int}(),
                               Vector{SubWorker{T,A,S}}(nworkers()),
                               Vector{Work}(nworkers()),
                               RemoteChannel(() -> DecisionChannel(Dict{Int,A}())),
                               RemoteChannel(() -> Channel{QCut{T}}(4*nworkers()*n)),
                               ξ₀_,
                               A(),
                               A(),
                               A(fill(-Inf,n)),
                               Vector{SparseHyperPlane{T}}(),
                               A(),
                               DRegularizedParameters{T}(;kw...),
                               ProgressThresh(1.0, "Distributed RD L-Shaped Gap "))
        lshaped.progress.thresh = lshaped.parameters.τ
        push!(lshaped.subobjectives,zeros(n))
        push!(lshaped.finished,0)
        push!(lshaped.Q_history,Inf)
        push!(lshaped.Q̃_history,Inf)
        push!(lshaped.θ_history,-Inf)

        init!(lshaped,subsolver)

        return lshaped
    end
end
DRegularized(model::JuMP.Model,mastersolver::AbstractMathProgSolver,subsolver::AbstractMathProgSolver; kw...) = DRegularized(model,rand(model.numCols),mastersolver,subsolver; kw...)

function (lshaped::DRegularized)()
    # Reset timer
    lshaped.progress.tfirst = lshaped.progress.tlast = time()
    # Start workers
    active_workers = init_workers!(lshaped)
    # Start procedure
    while true
        status = iterate!(lshaped)
        if status != :Valid
            close_workers!(lshaped,active_workers)
            return status
        end
    end
end


text:
using Images, ImageMagick
function fourier_op_out(N::Int64, xmax::Float64, output_file::String)
    n = [2*xmax * i / N for i = 1:N]
    k = n'
    #n = [-xmax + 2*xmax*i/N for i =0:N-1]
    #k = append!(n[1:div(N,2)],n[div(N,2):N-1])'
    transform_matrix = (exp.(-2im*pi*n*k/N).+(1.0+1im))/2
    a = Array{RGB{N0f8},2}(undef, N, N)
    for i = 1:N
        for j = 1:N
            a[i,j] = RGB{N0f8}(real(transform_matrix[i,j]), 0, imag(transform_matrix[i,j]))
        end
    end

    save(output_file, a)

end

function DFT_out(N::Int64, output_file::String)
    n = [i for i = 0:N-1]
    k = n'
    transform_matrix = (exp.(-2im*pi*n*k/N).+(1.0+1im))/2
    a = Array{RGB{N0f8},2}(undef, N, N)
    for i = 1:N
        for j = 1:N
            a[i,j] = RGB{N0f8}(real(transform_matrix[i,j]), 0, imag(transform_matrix[i,j]))
        end
    end

    save(output_file, a)

end


text:
@testset verbose = true "Smoketests" begin

    import DataFrames
    import Random: rand,randstring

##################################################
# Read the data
##################################################

    thisfile = joinpath(datadir, "GermanMotorPremiums", "data1Small.csv")
    @test isfile(thisfile)
    @time df_tmp = CSV.read(thisfile, DataFrame,strict=true, types=eltypesData1,  pool=true, copycols=true);
    df_tmp_orig = deepcopy(df_tmp)

# example with only char columns
    selected_explanatory_vars_char_only = ["ART_DES_WOHNEIGENTUM","FAMILIENSTAND","GESCHLECHT","FINANZIERUNGSART","STADT","KENNZEICHEN"]
    dtmtableMini, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp[1:100,:], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars_char_only);
    sett.minWeight = -.2
    @info("DTM:Testing build tree (dtm mini; only character independent variables)")
    strs, resm = dtm(dtmtableMini, sett)
    @test true

    settB = deepcopy(sett)
    settB.model_type = "boosted_tree"
    @info("DTM:Testing boosted tree (dtm mini; only character independent variables)")
    strs, resm = dtm(dtmtableMini, settB)
    @test true

# example with only numeric columns
    selected_explanatory_vars_num_only = ["GEBURTSDATUM","VORSCHAEDEN_ANZAHL","AUTOMOBILCLUB_MITGLIED_SEIT"]
    dtmtableMini, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp[1:100,:], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars_num_only);
    sett.minWeight = -.2
    @info("DTM:Testing build tree (dtm mini; only numerical independent variables)")
    strs, resm = dtm(dtmtableMini, sett)
    @test true

    settB = deepcopy(sett)
    settB.model_type = "boosted_tree"
    @info("DTM:Testing boosted tree (dtm mini; only numerical independent variables)")
    strs, resm = dtm(dtmtableMini, settB)

# other tests:
    selected_explanatory_vars = ["PLZ_WOHNORT","ART_DES_WOHNEIGENTUM","GEBURTSDATUM","FAMILIENSTAND","NATIONALITAET","GESCHLECHT","FINANZIERUNGSART","STADT","KENNZEICHEN"]

# keep 10 largest PLZ only
    vorig = deepcopy(df_tmp[!,:PLZ_WOHNORT])
    counts, freqs, vals, keep = getCounts(vorig, threshold=10)
    vnew = df_tmp[!,:PLZ_WOHNORT]
    mapToOther!(vnew, keep, 9999999)

##################################################
# Run Mini Tree
##################################################
    dtmtableMini, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp[1:100,:], treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars);
    sett.minWeight = -.2
    @info("DTM:Testing build tree (dtm mini)")
    strs, resm = dtm(dtmtableMini, sett)

# boosting
    settB = deepcopy(sett)
    settB.model_type = "boosted_tree"
    @info("DTM:Testing boosted tree (dtm mini)")
    strs, resm = dtm(dtmtableMini, settB)

############################################
# graphviz
############################################

    @testset verbose = true "Graphviz Smoketests" begin
        @show graphvizexe
        @show graphvizsetup_is_working
        if graphvizsetup_is_working
            #do graphviz tests
                dtmtableMini, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp[1:100,:], treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars,print_details=false);
                sett.write_dot_graph = true
                sett.graphvizexecutable = graphvizexe
                sett.minWeight = -.2
                strs, resm = dtm(dtmtableMini, sett);
                @test true
                pdf_files_produced = filter(x->endswith(lowercase(x),".pdf"),strs)
                @test size(pdf_files_produced,1)==1

            # boosting
                settB = deepcopy(sett)
                settB.model_type = "boosted_tree"
                @info("DTM:Testing boosted tree (dtm mini)")
                strs, resm = dtm(dtmtableMini, settB)
                @test true
                pdf_files_produced = filter(x->endswith(lowercase(x),".pdf"),strs)
                @test size(pdf_files_produced,1)==2

        else
            #GraphViz Exe not found
            @test false
            @warn("""Graphvizexecutable not found. Please specify the path to dot.exe in ENV["graphvizdot"]""")
        end
    end #graphviz smoketests

# selected_explanatory_vars=[    "VORSCHAEDEN_ANZAHL",    "MALLORCA_POLICE",	"SCHUTZBRIEF_INKL",	"FREIE_WERKSTATTWAHL",	"AUTOMOBILCLUB_MITGLIED_SEIT",	"BAHNCARD",	"ZAHLUNGSWEISE",	"JAHRESKARTE_OEPNV",	"MOTORRAD_BESITZER",	"AUTOMOBILCLUB",	"SFKLASSE_VOLLKASKO",	"SFKLASSE_HAFTPFLICHT",	"STELLPLATZ_ABSCHLIESSBAR",	"NAECHTLICHER_STELLPLATZ",	"NUTZUNGSWEISE",	"JAEHRLICHE_FAHRLEISTUNG",	"TSN",	"ERSTZULASSUNG",	"HSN",	"FINANZIERUNGSART",	"ZULASSUNG_AUF_VERSICHERUNGSNEHM",	"STADT",	"KENNZEICHEN",	"PLZ_DES_HALTER",	"SELBSTGENUTZTES_WOHNEIGENTUM",	"ART_DES_WOHNEIGENTUM",	"GEBURTSDATUM",	"FAMILIENSTAND",	"NATIONALITAET",	"GESCHLECHT",	"FUEHRERSCHEIN_ERWORBEN_AM",	"VORSCHAEDEN0_typeKH",	"VORSCHAEDEN0_typetk",	"VORSCHAEDEN0_month",	"VORSCHAEDEN0_year",	"VORSCHAEDEN1_typetk",	"VORSCHAEDEN1_month",	"VORSCHAEDEN1_year",	"VORSCHAEDEN2_typevk",	"VORSCHAEDEN2_month",	"VORSCHAEDEN2_year",	"adacid",	"name",	"marke",	"modell",	"preis",	"getriebeart",	"antriebsart",	"Fahrzeugklasse",	"co2klasse",	"kw",	"ps",	"tueranzahl",	"Motorart",	"Kraftstoffart",	"Motorbauart",	"Schadstoffklasse",	"Karosserie",	"Sitzanzahl",	"typklasseh_num",	"typklassetk_num",	"typklassevk_num",	"hubraum2",	"drehmoment2",	"breite2",	"radstand2",	"laenge2",	"hoehe2",	"leergewicht2",	"gesamtgewicht2",	"zuladung2",	"kofferraumvolumen_num",	"hoechstgeschwindigkeit2",	"verbrauchgesamt2",	"verbrauchausserorts2",	"verbrauchinnerorts2",	"beschleunigung2",	"tank2",	"kfzsteuer2",	"anzahlgaenge2",	"anzahlzylinder2",	"co2_wert",	"modellstart_y"]
##################################################
# run tree
##################################################
    dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars,print_details=false);
    @info("DTM:Testing build tree")
    strs, resm = dtm(dtmtable, sett)
    @test typeof(resm.modelstats) == DataFrame
    @test 1 == 1

# predict
    DecisionTrees.predict(resm, dtmtable.features)

##################################################
# run boosting
##################################################

    sett.iterations = 10
    sett.model_type = "boosted_tree"
    @info("DTM:Testing boosted tree. Iterations=$(sett.iterations)")
    strs, resm2 = dtm(dtmtable, sett)
    @test typeof(resm2.modelstats) == DataFrame
    @test 1 == 1

# reduce iterations for the remaining tests
    sett.iterations = 3

    sett.statsByVariables = Int[1,2]
    @info("DTM:Testing statsByVariables")
    strs, resm2 = dtm(dtmtable, sett)
    sett.statsByVariables = Int[]

    @info("DTM:Testing different splitting criterias...")

# predict
DecisionTrees.predict(resm2, dtmtable.features)

# try different splitting criteria
    for splitCrit in ["difference","poissondeviance","gammadeviance","mse","sse","msepointwise"], modelTYPE in ["build_tree","boosted_tree"]
        updateSettingsMod!(sett, crit=splitCrit, model_type=modelTYPE)
        try
            strs, resmT = dtm(dtmtable, sett)
            @test typeof(resmT.modelstats) == DataFrame
        catch thisErr
            @show stacktrace()
            @show thisErr
            @test "This one failed->" == string("crit=", splitCrit, "; model_type=", modelTYPE)
        end
    end

    # try different splitting criteria with randomw > 0
    for splitCrit in ["difference","poissondeviance","gammadeviance"], modelTYPE in ["build_tree","boosted_tree"]
        updateSettingsMod!(sett, crit=splitCrit, model_type=modelTYPE)
        updateSettingsMod!(sett, randomw = 0.05)
        try
            strs, resmT = dtm(dtmtable, sett)
            @test typeof(resmT.modelstats) == DataFrame
        catch thisErr
            @show stacktrace()
            @show thisErr
            @test "This one failed->" == string("crit=", splitCrit, "; model_type=", modelTYPE, "; randomw=", rw)
        end
    end

    updateSettingsMod!(sett, crit="difference")
    @info("DTM:Finished testing different splitting criterias.")

##################################################
# run bagging
##################################################
# tbd need to update bagging to work under j0.7 version of Decisiontrees
    sett.model_type = "bagged_tree"
# @info("DTM:testing bagging")
# dtm(dtmtable,sett)
# @test 1==1

##################################################
# check case where we have more than 255 splitting points
##################################################

    counts, freqs, vals, keep = getCounts(vorig, threshold=300)
    vnew = deepcopy(vorig)
    mapToOther!(vnew, keep, 9999999)
    df_tmp[!,:PLZ_WOHNORT] = vnew

    dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars,print_details=false);

    sett.maxSplittingPoints = 300
    dtmtable = prep_data_from_df(df_prepped, sett, joinpath(mktempdir(), "dtmsmoketest.csv"),print_details=false);
    @test eltype(dtmtable.features[!,:PLZ_WOHNORT].refs) == UInt16
    @info("DTM:testing maxSplittingPoints=$(sett.maxSplittingPoints)")
    strs, resm2b = dtm(dtmtable, sett)
    @test typeof(resm2b.modelstats) == DataFrame

##################################################
# More Tests "for coverage" purposes
##################################################

##################################################
# Multirun Boosting
##################################################
    settV = createGridSearchSettings(sett,
    minWeight=[-0.1,-0.2]
    ,learningRate=[0.1,0.05],iterations=[2])
    settVBoosting = deepcopy(settV)
    @info("DTM:testing gridsearch for boosting")
    a, b, allmodels = dtm(dtmtable, settV)
    @test typeof(allmodels[1].modelstats) == DataFrame

##################################################
# Multirun buildtree
##################################################
    sett.model_type = "build_tree"
    settV = createGridSearchSettings(sett,
    minWeight=[-0.1,-0.2]
    ,learningRate=[0.1,0.05],iterations=[2])
    @info("DTM:testing gridsearch/multirun for build_tree")
    a, b, allmodels = dtm(dtmtable, settV)
# @test typeof(allmodels[1].modelstats)==DataFrame

##################################################
# Cross validation
##################################################
    sett.model_type = "boosted_tree"
    cvsampler = CVOptions(-3, 0.0, true)
    @info("DTM:testing cvsampler")
    statsdf, settsdf, cvModels = dtm(dtmtable, sett, cvsampler)

    cvsampler = CVOptions(3, 0.65, true)
    @info("DTM:testing cvsampler2")
    statsdf, settsdf, cvModels = dtm(dtmtable, sett, cvsampler)

    @warn("We can possibly enable the following snippet now (it should run fine)")
    @warn("CV sampler tests were breaking the other tests (possibly because dtmtable is modified which should not happen!")
    if true # false
        @info("DTM:testing cvsampler3")
        cvsampler = CVOptions(3, 0.5, false)
        statsdf, settsdf, cvModels = dtm(dtmtable, sett, cvsampler)
    else
        @test_broken 1 == 2
    end

    @info("DTM:testing different performance measures")
    cvsampler = CVOptions(-3, 0.0, true)
    for thisPerfMeasure in DecisionTrees.global_statsperiter_header
        updateSettingsMod!(sett, performanceMeasure=thisPerfMeasure)
        statsdf = DataFrame()
        statsdf, settsdf, cvModels = dtm(dtmtableMini, sett, cvsampler)
        @test size(statsdf, 1) > 0
    end

# set some default settings
    updateSettingsMod!(sett,
iterations=3,
model_type="boosted_tree",
nScores="1000",
writeStatistics="true",
writeSasCode="true",
writeIterationMatrix="true",
writeResult="true",
writeCsharpCode="true",
writeVbaCode="true",
saveJLDFile="true",
saveResultAsJLDFile="false",

showProgressBar_time="true",
prroduceEstAndLeafMatrices="true",
write_dot_graph="false",
calculateGini="true",

calculatePoissonError="true",
performanceMeasure="Average Poisson Error Val"
)
    sett.model_type = "build_tree"
    @info("DTM:testing a build_tree model")
    strs, resm3 = dtm(dtmtable, sett)
    sett.model_type = "boosted_tree"
    @info("DTM:testing a boosted_tree")
    strs, resm3 = dtm(dtmtable, sett)
    @test typeof(resm3.modelstats) == DataFrame
# test if files exist
    for fi in strs
        @test isfile(fi)
    end

# more smoketests
    sett.subsampling_prop = 0.7
    sett.subsampling_features_prop = 0.7
    @info("DTM:testing subsampling_features_prop and subsampling_prop")
    strs, resm4 = dtm(dtmtable, sett)
    @test typeof(resm4.modelstats) == DataFrame

    sett.boolRandomizeOnlySplitAtTopNode = false
    sett.randomw = 0.02
    @info("DTM:testing randomw=$(sett.randomw)")
    strs, resm5 = dtm(dtmtable, sett)
    @test typeof(resm5.modelstats) == DataFrame

    sett.spawnsmaller = false
    sett.randomw = 0.0
    @info("DTM:testing spawnsmaller=$(sett.spawnsmaller) randomw=$(sett.randomw)")
    strs, resm5b = dtm(dtmtable, sett)
    @test typeof(resm5b.modelstats) == DataFrame

# minWeight too big
    oldminw = sett.minWeight
    sett.minWeight = -0.51
    @info("DTM:testing minWeight=$(sett.minWeight)")
    strs, resm6 = dtm(dtmtable, sett)
    @test typeof(resm6.modelstats) == DataFrame

# for coverage
# empty keycol
    @info("DTM:testing errors/warnings of prepare_dataframe_for_dtm")
    prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], keycol="", numcol="LOSS20HALF", independent_vars=selected_explanatory_vars);
    df_tmp2 = deepcopy(df_tmp_orig)
    df_tmp2[!,:PREMIUM66][1:20] .= -20.1
    dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars,print_details=false);
    df_tmp2[!,:PREMIUM66][1:20] .= 0.0
    dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars,print_details=false);
    df_tmp2[!,:LOSS20HALF][1:20] .= 0.0
    dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars,print_details=false);
    df_tmp2[!,:LOSS20HALF][1:20] .= -20
    dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars,print_details=false);

    sett.minWeight = oldminw
# example with weights.==1 and denominator.==1
    dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], numcol="LOSS20HALF", independent_vars=selected_explanatory_vars,print_details=false);

# multithreaded runs
    if false  # currently disabled
        @info("DTM:testing multithreading")
        @test Distributed.nprocs() == 1 # we generally expect that tests are run on only 1 thread
        dtmtable, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp, treat_as_categorical_variable=["PLZ_WOHNORT"], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars,print_details=false);
        sett.model_type = "boosted_tree"
        if Distributed.nprocs() < 2
            Distributed.addprocs(2)
            Distributed.@everywhere using Decisiontrees
        end
        statsdf, settsdf, cvModels = dtm(dtmtable, sett, cvsampler)
    # same for multirun
        a, b, allmodels = dtm(dtmtable, settV)
        dtm(dtmtable, settVBoosting)
    # remove workers again
        Distributed.rmprocs(workers())
    end

# coverage test of absrel
    DecisionTrees.absrel(rand(3), rand(3))

end # testset


text:
@testset "basic interest rates" begin

    ## vector interest rate
    @testset "vector interest rate" begin
        i = InterestRate([0.05, 0.05, 0.05])

        @test disc(i,0) == 1.0
        @test disc(i,1) == 1 / 1.05
        @test disc(i,2) == 1 / 1.05 ^ 2
        @test disc(i,3) == 1 / 1.05 ^ 3
        @test disc(i,1,2) == 1 / 1.05
        @test rate(i,1) == 0.05
        @test rate(i,2) == 0.05

        @test disc.(i,1:3) == [1 / 1.05 ^ t for t in 1:3]
    end

    ## real interest rate
    @testset "constant interest rate" begin
        i = InterestRate(0.05)

        @test disc(i,0) == 1.0
        @test disc(i,1) == 1 / 1.05
        @test disc(i,2) == 1 / 1.05 ^ 2
        @test disc(i,3) == 1 / 1.05 ^ 3
        @test disc(i,1,2) == 1 / 1.05
        @test rate(i,1) == 0.05
        @test rate(i,2) == 0.05

        @test disc.(i,1:3) == [1 / 1.05 ^ t for t in 1:3]
    end
end

@testset "DiscountFactor Iterator" begin

    target = [1 / 1.05 ^ t for t in 0:4 ]

    @testset "contsant" begin
        df = DiscountFactor(InterestRate(0.05),1)

        ds = Iterators.take(df,5) |> collect

        @test all(ds .≈ target)

        @test length(Iterators.take(df,100)) == 100

        # can't collect an infinite series
        @test_throws MethodError collect(df)

        @test df == InterestRate(0.05)(1)
    end

    @testset "vector" begin
    df = DiscountFactor(InterestRate(repeat([0.05],4)),1)

        ds = Iterators.take(df,5) |> collect
        @test all(ds .≈ target)
        @test all(collect(df) .≈ target)
        @test length(df) == 5

        @test_broken InterestRate([0.05,0.05,0.05,0.05])(1) == df #https://stackoverflow.com/questions/62336686/struct-equality-with-arrays
    end

end

@testset "passthrough life contingency" begin
    ins = LifeContingency(
        SingleLife(
            mort = UltimateMortality([0.5]),
            issue_age = 0
        ),
        InterestRate(0.05)
    )

    @test disc(ins,1) ≈ 1/1.05
    @test disc(ins,1,2) ≈ 1/1.05
end


text:
using DelimitedFiles

player_inv = Dict(
    "X" => "O",
    "O" => "X"
)

function choose_player_input()

    input = ""
    while true
        println("Enter which piece you want to play as (X or O):")
        input = readline()
        input = uppercase(input)
        println("Input chosen as $input")
        if (input == "X") | (input == "O")
            break
        else
            println("Please choose either X or O")
        end
    end
    return input
end

function have_we_won(board, ai)

    havewon = false
    winner = checkwinner(board)
    havewon = ai == winner ? true : false

    return havewon

end

function whosplaying(depth, scoreprofile)
    # since we start with ai. every odd depth is ai's turn
    # i.e. the turn of teh score profile with score +1
    if depth % 2 == 1
        player = [k for (k,v) in scoreprofile if v == 1][1]
    else
        player = [k for (k,v) in scoreprofile if v == -1][1]
    end

    return player

end

function minimax(board, scoreprofile; maximizer, depth = 1)
    #==
    score profile is a dictionary which shows whether X or O
    gets the +1 and -1
    ==#

    bestscore = nothing
    winner = checkwinner(board)
    player = whosplaying(depth, scoreprofile)
    if winner != ""
        bestscore = scoreprofile[winner]
    else
        possiblemoves = legal_moves(board)

        if maximizer
            bestscore = -10
            for move in possiblemoves
                if islegal(board[move...])
                    boardcopy = copy(board)
                    boardcopy[move...] = player
                    score = minimax(boardcopy, scoreprofile, maximizer = false, depth = depth + 1)
                    bestscore = max(score, bestscore)
                end

            end
        else
            bestscore = 10
            for move in possiblemoves
                if islegal(board[move...])
                    boardcopy = copy(board)
                    boardcopy[move...] = player
                    score = minimax(boardcopy, scoreprofile, maximizer = true, depth = depth + 1)
                    bestscore = min(score, bestscore)
                end

            end
        end

    end

    return bestscore
end

function ai_move!(board, ai)

    possiblemoves = legal_moves(board)
    scoreprofile = Dict(
        ai => 1,
        player_inv[ai] => -1,
        "Tie" => 0
    )
    bestscore = -10
    bestmove = (nothing, nothing)
    for move in possiblemoves

        boardcopy = copy(board)
        boardcopy[move...] = ai
        score = minimax(boardcopy, scoreprofile, maximizer = false, depth = 2)
        if score > bestscore
            bestscore = score
            bestmove = move
        end
    end

    board[bestmove...] = ai

end

function playgame_ai(nrow=3, ncol=3)

    board = create_board(nrow, ncol)
    nmaxturns = nrow * ncol
    gameoutcome = "Tie"

    human = choose_player_input()
    ai = player_inv[human]

    for i in 1:nmaxturns

        println("==========================")
        println("Turn number: $i")
        writedlm(stdout, board,  "\t|\t")

        player = whos_turn(board)

        if player == human
            play_legal_move!(board)
        elseif player == ai
            ai_move!(board, ai)
        else
            println("Something went wrong with choosing which player plays")
        end

        gameover = isgameover(board)
        if gameover
            winner = checkwinner(board)

            if winner == "Tie"
                gameoutcome = "Tie!"
            else
                gameoutcome = "$winner wins!"
            end

            break
        end

    end

    println("==========================")
    writedlm(stdout, board,  "\t|\t")

    return gameoutcome

end

max_stars_repo_name: HariharanJayashankar/tictactoe

max_issues_repo_name: HariharanJayashankar/tictactoe

max_forks_repo_name: HariharanJayashankar/tictactoe