text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    lept_version()::String

Return the version number of the leptonica library.
"""
function lept_version()::String
    local retval = ccall(
        (:getLeptonicaVersion, LEPTONICA),
        Ptr{UInt8},
        ()
    )
    local result = ""

    if retval != C_NULL
        result = unsafe_string(retval)
        lept_free(retval)
    end

    return result
end

# =================================================================================================
"""
    lept_free(
        ptr::Ptr
    )::Nothing

Call the leptonica library to free some memory it allocated.

__Arguments:__
| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | ptr  |         | The pointer to the memory allocated by leptonica.

__Details:__

This method is meant for internal usage.  At no point is memory allocated by leptonica passed to
clients of this library.
"""
function lept_free(
            ptr::Ptr
        )::Nothing
    ccall(
        (:lept_free, LEPTONICA),
        Cvoid,
        (Ptr{Cvoid},),
        ptr
    )
    nothing
end


text:
@doc """
    start_recording()

Start recording.

See also: [`stop_recording`](@ref)
"""
start_recording

@doc """
    stop_recording()

Stop recording.

See also: [`start_recording`](@ref)
"""
stop_recording

@doc """
    pause_recording()

Pause recording of active recorder.

See also: [`restart_recording`](@ref)
"""
pause_recording

@doc """
    restart_recording()

Restart recording of last active recorder.

See also: [`pause_recording`](@ref)
"""
restart_recording

let
    global recorder, start_recording, stop_recording, pause_recording, restart_recording
    _recorders::Dict{String, Base.Process}                 = Dict{String, Base.Process}()
    _active_recorder_id::String                            = ""
    recorder(id::String=DEFAULT_RECORDER_ID)::Base.Process = _recorders[id]
    active_recorder_id()::String                           = (if (_active_recorder_id=="") error("no recorder is active.") end; return _active_recorder_id)

    function start_recording(id::String=DEFAULT_RECORDER_ID)
        _recorders[id] = open(`$RECORDER_BACKEND $RECORDER_ARGS`) # NOTE: here could be started multiple recorders.
        _active_recorder_id = id
    end

    function stop_recording(id::String=DEFAULT_RECORDER_ID)
        close(_recorders[id])
    end

    pause_recording()   = (stop_recording(active_recorder_id()); return)
    restart_recording() = (start_recording(active_recorder_id()); return)
end

ext: jl


text:
#This file contains all the tests of the Utils Functions
using Base.Test
include("LoadUtils.jl")

@testset "Utility Functions Test" begin
    @test datetime2modifiedjulian(DateTime(2016,1,1,0,0,5,0)) == datetime2modifiedjulian(DateTime(2016,1,1,0,0,0,0),5)
end

##Examples
#x = [0 pi/4 pi/2 3*pi/4 pi 5*pi/4 6*pi/4 7*pi/4 2*pi]
#y = x
#z = y
#v = sin.(x)
#v2 = repmat(sin.(x),9,1) + repmat(sin.(y)',1,9)

#vq = lininterpn(x,v,1)
#vqq = lininterp1(x,v,1)
#vq2 = lininterpn(x,y,v2,1,1)
#vqq2 = lininterp2(x,y,v2,1,1)


text:
module Intermediate

#
#   This module provides some intial calculations for things like numbers and
#   ages of children, number of carers and so on, for benefit units and households.
#   It is intended to simplify subsequent modelling of benefits, etc..
#   and to reduce the number of repeated calculations. It also holds code for some age and
#   child count calculations that are also used stand-alone elsewhere.
#
#   FIXME some of the names should be globally changed.
#

using Base: Bool, String
using Dates: TimeType, Date, now, Year

using ScottishTaxBenefitModel
using .Definitions

using .ModelHousehold:
    BenefitUnit,
    Household,
    Person,
    between_ages,
    count,
    default_bu_allocation,
    empl_status_in,
    ge_age,
    get_benefit_units,
    get_head,
    has_carer_member,
    has_children,
    has_disabled_member,
    has_income,
    is_lone_parent,
    is_severe_disability,
    is_single,
    le_age,
    num_adults,
    num_carers,
    num_people,
    pers_is_carer,
    pers_is_disabled,
    search

using .STBParameters:
    AgeLimits,
    HoursLimits,
    ChildLimits,
    reached_state_pension_age

using .Utils:
    has_non_z,
    haskeys,
    mult,
    to_md_table

using .STBIncomes

export
    MTIntermediate,
    apply_2_child_policy,
    born_before,
    has_limited_capactity_for_work_activity,
    has_limited_capactity_for_work,
    is_working_hours,
    make_intermediate,
    make_recipient,
    num_born_before,
    to_string,
    working_disabled,
    HHIntermed

"""
examples:

2 children born before start_date, 1 after
 => allowable = 2
3  children born before start_date, 1 after
 => allowable = 3
1  children born before start_date, 2 after
 => allowable = 2
0  children born before start_date, 2 after
 => allowable = 2
 @return number of children allowed
"""
function apply_2_child_policy(
    bu      :: BenefitUnit,
    limits  :: ChildLimits
    ;
    model_run_date :: TimeType = now() ) :: Integer
    before_children = 0
    after_children = 0
    for pid in bu.children
        ch = bu.people[pid]
        if born_before( ch.age, limits.policy_start, model_run_date )
            before_children += 1
        else
            after_children += 1
        end
    end
    allowable = before_children + min( max(limits.max_children-before_children,0), after_children )
end

function born_before( age :: Integer,
    start_date     :: TimeType = Date( 2017, 4, 6 ), # 6th April 2017
    model_run_date :: TimeType = now() )
    bdate = model_run_date - Year(age)
    return bdate < start_date
end

"""
Very loose implementation of CPAG 2020/1 ch 45
"""
function has_limited_capactity_for_work( pers :: Person ) :: Bool
    # based on cpag 2020/1 ch45 pp 1100-
    l =  pers.employment_status in
        [   Retired,
            Permanently_sick_or_disabled,
            Temporarily_sick_or_injured,
            Other_Inactive] &&  # FIXME - retired/temp sick/looking after family ??
        pers_is_disabled( pers )
    l = l || (pers.esa_type != no_jsa)
    return l
end

"""
Very loose implementation of CPAG 2020/1 ch 45. Until we have
something better, maybe regression-based.
`Work Activity` here means attending interviews, assessments and the like.
"""
function has_limited_capactity_for_work_activity( pers :: Person ) :: Bool
    l = false
    if has_limited_capactity_for_work( pers )
        has_high_esa = false
        if has_non_z( pers.income, employment_and_support_allowance )
            esa = pers.income[employment_and_support_allowance]
            if( pers.age < 25 )
                if esa > 80 # crudely including support component
                    has_high_esa = true
                end
            else
                if esa > 100 # crudely including support component
                    has_high_esa = true
                end
            end
        end # crude ESA check
        if is_severe_disability( pers ) || has_high_esa
            l = true
        end
    end
    return l
end

function num_born_before(
    bu             :: BenefitUnit,
    start_date     :: TimeType = Date( 2017, 4, 6 ), # 6th April 2017
    model_run_date :: TimeType = now()) :: Integer
    nb = 0
    for pid in bu.children
        ch = bu.people[pid]
        if born_before( ch.age, start_date, model_run_date )
            nb += 1
        end
    end
    return nb
end

function is_working_hours( pers :: Person, hours... ) :: Bool
    if length(hours) == 1
        return (pers.usual_hours_worked >= hours[1])
    elseif length(hours) == 2
        return (pers.usual_hours_worked >= hours[1]) &&
               (pers.usual_hours_worked <= hours[2])
    end
end

function is_working_lt_hours( pers :: Person, hours :: Real ) :: Bool
    return pers.usual_hours_worked < hours
end

"""
See CPAG ch 61 p 1426 and appendix 5
"""
function working_disabled( pers::Person, hrs :: HoursLimits ) :: Bool
    if pers.usual_hours_worked >= hrs.lower || pers.employment_status in [Full_time_Employee, Full_time_Self_Employed]
        if pers.registered_blind || pers.registered_partially_sighted || pers.registered_deaf
            return true
        end
        for (dis, t ) in pers.disabilities
            return true
        end
        if haskeys( pers.income,
            [
                Incapacity_Benefit,
                Severe_Disability_Allowance,
                Employment_and_Support_Allowance ])
            return true
        end
    end
    return false
end

#
# FIXME some names here need clarified
#
mutable struct MTIntermediate

    benefit_unit_number :: Int
    num_people :: Int
    age_youngest_adult :: Int
    age_oldest_adult :: Int
    age_youngest_child :: Int
    age_oldest_child :: Int
    num_adults :: Int
    someone_pension_age  :: Bool
    someone_pension_age_2016 :: Bool
    all_pension_age :: Bool
    someone_working_ft  :: Bool
    #
    someone_working_ft_and_25_plus :: Bool

    num_not_working :: Int
    num_working_ft :: Int
    num_working_pt :: Int
    num_working_24_plus :: Int
    num_working_16_or_less :: Int
    total_hours_worked :: Int
    someone_is_carer :: Bool
    num_carers :: Int

    is_sparent  :: Bool
    is_sing  :: Bool # FIXME RENAME: is_single
    is_disabled :: Bool

    num_disabled_adults :: Int
    num_disabled_children :: Int
    num_severely_disabled_adults :: Int
    num_severely_disabled_children :: Int

    num_job_seekers :: Int

    num_children :: Int
    num_allowed_children :: Int
    num_children_born_before :: Int
    ge_16_u_pension_age  :: Bool
    limited_capacity_for_work  :: Bool
    has_children  :: Bool
    economically_active :: Bool
    working_disabled :: Bool
    num_benefit_units :: Int
end

struct HHIntermed
    hhint :: MTIntermediate
    buint ::  Vector{MTIntermediate}
end

function to_string( it :: MTIntermediate )::String
    s = to_md_table( it, depth=2 )
    return s
end

function to_string( hh :: HHIntermed ) :: String
    s = "## Household\n"
    s *= to_string( hh.hhint )
    for i in eachindex( hh.buint )
        s *= "### Bu $i"
        s *= to_string( hh.buint[i])
    end
    return s
end

"""
Assign CB, etc. to 1st adult female in a benefit unit, everything else to head of bu
"""
function make_recipient( bu :: BenefitUnit, which :: Incomes ) :: BigInt
    ## FIXME assert check here for non bu level benefits e.g. PIPs
    #
    # FIXME maybe also check relationships to children?
    if which in [CHILD_TAX_CREDIT, CHILD_BENEFIT, SCOTTISH_CHILD_PAYMENT ]
        # first_adult_female
        for pid in bu.adults
            if bu.people[pid].sex == Female
                return pid
            end
        end
    end
    ## FIXME expand this to check for if head is retired for WTC/UC
    # see the thing in .UniversalCredit.jl
    return get_head( bu ).pid
end

function aggregate!( sum :: MTIntermediate, add :: MTIntermediate )
    # benefit_unit_number :: Int
    sum.num_people += add.num_people
    sum.age_youngest_adult = min( sum.age_youngest_adult, add.age_youngest_adult )
    sum.age_oldest_adult = max( sum.age_oldest_adult, add.age_oldest_adult )
    sum.age_youngest_child = min( sum.age_youngest_child, add.age_youngest_child )
    sum.age_oldest_child = max( sum.age_oldest_child, add.age_oldest_child )
    sum.num_adults += add.num_adults
    sum.someone_pension_age = add.someone_pension_age || add.someone_pension_age
    sum.someone_pension_age_2016 = sum.someone_pension_age_2016 || add.someone_pension_age_2016
    sum.all_pension_age = sum.all_pension_age && add.all_pension_age
    sum.someone_working_ft  =  sum.someone_working_ft  || add.someone_working_ft # someone FIXME
    sum.someone_working_ft_and_25_plus =  sum.someone_working_ft_and_25_plus  || add.someone_working_ft_and_25_plus # someone FIXME

    sum.num_working_pt += add.num_working_pt
    sum.num_working_ft += add.num_working_ft
    sum.num_not_working += add.num_not_working

    sum.num_working_24_plus += add.num_working_24_plus
    sum.num_working_16_or_less += add.num_working_16_or_less
    sum.total_hours_worked += add.total_hours_worked
    sum.someone_is_carer = sum.someone_is_carer || add.someone_is_carer # rename someone FIXME
    sum.num_carers += add.num_carers

    sum.is_sparent = sum.is_sparent || add.is_sparent
    sum.is_sing = false # since we're adding a 2nd bu
    sum.is_disabled = sum.is_disabled || add.is_disabled

    sum.num_disabled_adults += add.num_disabled_adults
    sum.num_disabled_children += add.num_disabled_children
    sum.num_severely_disabled_adults += add.num_severely_disabled_adults
    sum.num_severely_disabled_children += add.num_severely_disabled_children

    sum.num_job_seekers += add.num_job_seekers

    sum.num_children += add.num_children
    sum.num_allowed_children += add.num_allowed_children
    sum.num_children_born_before += add.num_children_born_before
    sum.ge_16_u_pension_age = sum.ge_16_u_pension_age || add.ge_16_u_pension_age # SOMEONE: FIXME RENAME
    sum.limited_capacity_for_work = sum.limited_capacity_for_work || add.limited_capacity_for_work
    sum.has_children = sum.has_children || add.has_children
    sum.economically_active = sum.economically_active || add.economically_active
    sum.working_disabled = sum.working_disabled || add.working_disabled
end

function make_intermediate(
    buno :: Int,
    bu   :: BenefitUnit,
    hrs  :: HoursLimits,
    age_limits :: AgeLimits,
    child_limits :: ChildLimits,
    num_benefit_units :: Int ) :: MTIntermediate
    # {RT} where RT
    age_youngest_adult :: Int = 9999
    age_oldest_adult :: Int = -9999
    age_youngest_child :: Int = 9999
    age_oldest_child :: Int = -9999

    is_working_disabled :: Bool = false
    num_adlts :: Int = num_adults( bu )
    # check both & die as we need to think about counts again
    # if we're going back in time to when these were unequal
    #
    num_pens_age :: Int = 0
    ge_16_u_pension_age  :: Bool = false
    someone_pension_age_2016 :: Bool = false
    limited_capacity_for_work = false
    num_working_16_or_less  = 0
    # :: Int = count( bu, is_working_lt_hours, hrs.lower )

    for pid in bu.adults
        pers = bu.people[pid]
        if has_limited_capactity_for_work( pers )
            limited_capacity_for_work = true
        end
        if pers.usual_hours_worked < hrs.lower
            num_working_16_or_less += 1
        end
        if reached_state_pension_age(
            age_limits,
            pers.age,
            pers.sex )
            num_pens_age += 1
        else
            ge_16_u_pension_age = true
        end
        if reached_state_pension_age(
            age_limits,
            pers.age,
            pers.sex,
            age_limits.savings_credit_to_new_state_pension )
            someone_pension_age_2016 = true
        end

    end
    someone_pension_age  :: Bool = num_pens_age > 0
    all_pension_age :: Bool = num_adlts == num_pens_age
    num_working_ft :: Int = count( bu, is_working_hours, hrs.higher )
    num_working_pt :: Int = count( bu, is_working_hours, hrs.lower, hrs.higher-0.000001 ) # eps(FloatType)
    num_not_working :: Int = count( bu, empl_status_in,
       Unemployed,
       Retired,
       Student,
       Looking_after_family_or_home,
       Permanently_sick_or_disabled,
       Temporarily_sick_or_injured,
       Other_Inactive
        )

    num_job_seekers = 0
    num_working_24_plus :: Int = count( bu, is_working_hours, hrs.med )
    someone_working_ft  :: Bool = num_working_ft > 0
    someone_working_ft_and_25_plus :: Bool = false

    total_hours_worked :: Int = 0
    num_carrs :: Int = num_carers( bu )
    someone_is_carer :: Bool = has_carer_member( bu )
    is_sparent  :: Bool = is_lone_parent( bu )
    is_sing  :: Bool = is_single( bu )

    has_children  :: Bool = ModelHousehold.has_children( bu )
    economically_active = search( bu, empl_status_in,
        Full_time_Employee,
        Part_time_Employee,
        Full_time_Self_Employed,
        Part_time_Self_Employed,
        Unemployed,
        Temporarily_sick_or_injured )

    # can't think of a simple way of doing the rest with searches..

    is_disabled = has_disabled_member( bu )
    num_disabled_adults = 0
    num_severely_disabled_adults = 0
    for pid in bu.adults
        pers = bu.people[pid]

        if( pers.age >= 25 ) && ( pers.usual_hours_worked >= hrs.higher )
            someone_working_ft_and_25_plus = true
        end

        if (pers.employment_status == Unemployed) || ( pers.jsa_type != no_jsa )
            num_job_seekers += 1
        end

        if pers.age > age_oldest_adult
            age_oldest_adult = pers.age
        end
        if pers.age < age_youngest_adult
            age_youngest_adult = pers.age
        end
        total_hours_worked += round(pers.usual_hours_worked)
        if working_disabled( pers, hrs )
            is_working_disabled = true
        end
        if pers_is_disabled( pers )
            num_disabled_adults += 1
            if is_severe_disability( pers )
                num_severely_disabled_adults += 1
            end
        end
    end
    @assert 120 >= age_oldest_adult >= age_youngest_adult >= 16
    num_children = size( bu.children )[1] # count( bu, le_age, 16 )
    num_children_born_before = num_born_before( bu ) # fixme parameterise
    num_disabled_children = 0
    num_severely_disabled_children :: Int = 0
    for pid in bu.children
        pers = bu.people[pid]
        if pers.age > age_oldest_child
            age_oldest_child = pers.age
        end
        if pers.age < age_youngest_child
            age_youngest_child = pers.age
        end
        if pers_is_disabled( pers )
            num_disabled_children += 1
            if is_severe_disability( pers )
                num_severely_disabled_children += 1
            end
        end
    end

    ## fixme parameterise this
    num_allowed_children :: Int = apply_2_child_policy( bu, child_limits )
    @assert (!has_children)||(19 >= age_oldest_child >= age_youngest_child >= 0)

    return MTIntermediate(
        buno,
        num_people( bu ),
        age_youngest_adult,
        age_oldest_adult,
        age_youngest_child,
        age_oldest_child,
        num_adlts,
        someone_pension_age,
        someone_pension_age_2016,
        all_pension_age,
        someone_working_ft,

        someone_working_ft_and_25_plus,

        num_not_working,
        num_working_ft,
        num_working_pt,
        num_working_24_plus,
        num_working_16_or_less,
        total_hours_worked,
        someone_is_carer,
        num_carrs,
        is_sparent,
        is_sing,
        is_disabled,
        num_disabled_adults,
        num_disabled_children,
        num_severely_disabled_adults,
        num_severely_disabled_children,

        num_job_seekers,

        num_children,
        num_allowed_children,
        num_children_born_before,
        ge_16_u_pension_age,
        limited_capacity_for_work,
        has_children,
        economically_active,
        is_working_disabled,
        num_benefit_units
    )
end

function make_intermediate(
    hh   :: Household,
    hrs  :: HoursLimits,
    age_limits :: AgeLimits,
    child_limits :: ChildLimits,
    allocator :: Function=default_bu_allocation ) :: HHIntermed

    bus = get_benefit_units( hh, allocator )
    n = size( bus )[1]
    buint = Vector{MTIntermediate}(undef,n)
    for buno in 1:n
        buint[buno] = make_intermediate(
            buno,
            bus[buno],
            hrs,
            age_limits,
            child_limits,
            n )
    end
    hhint = deepcopy( buint[1] )
    for buno in 2:n
        aggregate!( hhint, buint[buno])
    end
    return HHIntermed( hhint, buint )
end

end # module Intermediate

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using Test
using LogicCircuits
using LogicCircuits.IO: parse_lc_file, save_lines, parse_psdd_file, parse_sdd_file, parse_vtree_file

@testset "Circuit line parse-save-parse test" begin
    mktempdir() do tmp

        lines = parse_lc_file(zoo_lc_file("little_4var.circuit"))
        save_lines("$tmp/temp.circuit", lines)
        lines2 = parse_lc_file("$tmp/temp.circuit")
        @test length(lines) == length(lines2)

        lines = parse_psdd_file(zoo_psdd_file("little_4var.psdd"))
        save_lines("$tmp/temp.psdd", lines)
        lines2 = parse_psdd_file("$tmp/temp.psdd")
        @test length(lines) == length(lines2)

        lines = parse_sdd_file(zoo_sdd_file("random.sdd"))
        save_lines("$tmp/temp.sdd", lines)
        lines2 = parse_sdd_file("$tmp/temp.sdd")
        @test length(lines) == length(lines2)

        lines = parse_vtree_file(zoo_vtree_file("little_4var.vtree"))
        # TODO split up vtree saver into decompilation phase and file writing phase so we can test like above

    end
    @test true
end


text:
# SPDX-FileCopyrightText: 2021 Johannes Hampp <johannes.hampp@zeu.jlug.de>
#
# SPDX-License-Identifier: MIT

# Activate julia environment
using Pkg
Pkg.activate("./workflow/envs")
using GlobalEnergyGIS

region = snakemake.wildcards["region"]
ssp = snakemake.wildcards["ssp"]
rcp = snakemake.wildcards["rcp"]
ssp_year = parse(Int64, snakemake.wildcards["ssp_year"])
era_year = parse(Int64, snakemake.wildcards["era_year"])

ssp_rcp_scenario = replace("$ssp-$rcp", "." => "") # GEGIS expects RCP without "."

# Use GEGIS ML approach to generated synthetic demand
predictdemand(gisregion=region,
        sspscenario=ssp_rcp_scenario,
        sspyear=ssp_year,
        era_year=era_year)

# Make sure output directory exists
mkpath(dirname(snakemake.output[1]))

# Move the output file from the default output location to a better location
mv("./resources/gegis/output/SyntheticDemand_" * region
        * "_" * ssp_rcp_scenario
        * "-" * string(ssp_year)
        * "_" * string(era_year)
        * ".jld",
   snakemake.output[1])

max_stars_repo_name: euronion/synde

max_issues_repo_name: euronion/synde

max_forks_repo_name: euronion/synde


text:
function benchmark(bd::BenchmarkData, method::RegressionMethod, nfolds::Int=5)

    nrange = length(bd.n)
    results = zeros(nrange, 6)

    validation_results = Array{ValidationResults,2}(nrange, nfolds)

    for i = 1:nrange
        # Generate some synthetic data
        rd = getdata(bd, i)
        true_support = find(abs.(rd.w) .> 1e-6)

        folds = kfolds(shuffleobs((rd.X', rd.Y)), k = nfolds)
        j = 0
        for ((X_train, y_train), (X_test, y_test)) in folds
            j += 1
            # Validate the right method parameters
            X_train, X_test = X_train', X_test'
            v_results = validate_params!(X_train, y_train, bd.sparsity, method)
            validation_results[i, j] = v_results
            tic()
            indices, w = solve_problem(method, X_train, y_train, bd.sparsity)
            time = toq()
            a = accuracy(indices, true_support)
            f = falsepositive(indices, true_support)
            pred_train = predict_sparse(X_train, indices, w)
            train_R2 = isRsquared(pred_train, y_train)
            pred_test  = predict_sparse(X_test, indices, w)
            test_R2 = oosRsquared(pred_test, y_test, y_train)
            @show bd.n[i], a
            results[i, :] .+= [a, f, train_R2, test_R2, time, method.gamma]
        end
    end
    results / nfolds, validation_results
end


text:
module TwistedHessianEC

export EllipticCurve, EllipticCurveElem, TwiHesEC, TwiHesECElem, identity, isIdentity, normalized,
jInvariant, Isogeny, allPointsNaive

include("FieldArithmetic.jl")

using .FieldArithmetic
import Base: +,-,*,/,^,==,inv,show,zeros,gcd,gcdx,iszero,isone,zero,one,identity
import .FieldArithmetic: order

export Ring, RingElem, CRing, CRingElem, Field, FieldElem, Zmod, ZmodElem, PolyRing, PolyRingElem,
QuotientRing, QuotientRingElem, Fpsquare, order, squareroot, cuberoot, cuberootOfOne, allcuberoots,
allElements

abstract type EllipticCurve end
abstract type EllipticCurveElem end

struct TwiHesEC{T <: CRing, S <: CRingElem} <: EllipticCurve
    parent::T
    a::S
    d::S
end

struct TwiHesECElem{T <: CRing, S <: CRingElem} <: EllipticCurveElem
    X::S
    Y::S
    Z::S
    curve::TwiHesEC{T,S}
end

# Curve-Constructor
function TwiHesEC(R::T, a::S, d::S) where {T <: CRing, S <: CRingElem}
    if isSingular(R,a,d)
        error("Curve is singular!")
    end
    return TwiHesEC{T,S}(R, a, d)
end

# Element-Constructor
function (E::TwiHesEC{T,S})(X::S, Y::S, Z::S) where {T <: CRing, S <: CRingElem}
    if !validPoint(X, Y, Z, E)
        error("Point ([$(X)] : [$(Y)] : [$(Z)]) is not valid!!!")
    end
    return TwiHesECElem{T,S}(X, Y, Z, E)
end

###### Printing ######
function show(io::IO, E::TwiHesEC{T,S}) where {T <: CRing, S <: CRingElem}
    print(io, "Elliptic curve [$(E.a)]*X^3 + Y^3 + Z^3 = [$(E.d)]*X*Y*Z defined over $(E.parent)")
end

function show(io::IO, P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    print(io, "([$(P.X)] : [$(P.Y)] : [$(P.Z)])")
end

###### Helper-functions ######
function isSingular(E::T, a::S, d::S) where {T <: CRing, S <: CRingElem}
    return (iszero(a*(27*a - d^3)))
end

function validPoint(X::S, Y::S, Z::S, E::TwiHesEC{T,S}) where {T <: CRing, S <: CRingElem}
    return E.a*X^3 + Y^3 + Z^3 == E.d*X*Y*Z && !isZeroCoordinate(X,Y,Z)
end

function isZeroCoordinate(X::S, Y::S, Z::S) where {S <: CRingElem}
    F = X.ring
    return (X,Y,Z)==(zero(F),zero(F),zero(F))
end

function identity(E::TwiHesEC{T,S}) where {T <: CRing, S <: CRingElem}
    F = E.parent
    return E(zero(F),-one(F),one(F))
end

function isIdentity(P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    return P == identity(P.curve)
end

function sameCurve(P::TwiHesECElem{T,S}, Q::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    return P.curve == Q.curve
end

function normalized(P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    X, Y, Z = P.X, P.Y, P.Z
    E = P.curve
    F = E.parent
    if !iszero(Z)
        Zinv = inv(Z)
        return E(X*Zinv, Y*Zinv, one(F))
    else
        Xinv = inv(X)
        return E(one(F), Y*Xinv, zero(F))
    end
end

#Snail speed check...
function order(P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    i = 1
    Pi = P
    while !isIdentity(Pi)
        Pi = Pi + P
        i += 1
    end
    return i
end

###### Basic operators #######

function -(P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    E = P.curve
    return E(P.X, P.Z, P.Y)
end

function +(P::TwiHesECElem{T,S}, Q::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    if !(sameCurve(P,Q))
        error("Cannot add points from different curves!")
    end
    E = P.curve
    X1, Y1, Z1 = P.X, P.Y, P.Z
    X2, Y2, Z2 = Q.X, Q.Y, Q.Z
    X = X1^2*Y2*Z2 - X2^2*Y1*Z1
    Y = Z1^2*X2*Y2 - Z2^2*X1*Y1
    Z = Y1^2*X2*Z2 - Y2^2*X1*Z1
    if isZeroCoordinate(X,Y,Z)
        return fallbackAdd(P,Q)
    else
        return E(X, Y, Z)
    end
end

function fallbackAdd(P::TwiHesECElem{T,S}, Q::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    E = P.curve
    a = E.a
    X1, Y1, Z1 = P.X, P.Y, P.Z
    X2, Y2, Z2 = Q.X, Q.Y, Q.Z
    X = Z2^2*X1*Z1 - Y1^2*X2*Y2
    Y = Y2^2*Y1*Z1 - a*X1^2*X2*Z2
    Z = a*X2^2*X1*Y1 - Z1^2*Y2*Z2
    return E(X, Y, Z)
end

function -(P::TwiHesECElem{T,S}, Q::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    return P+(-Q)
end

#Regard E as a Z-module... This is again square and multiply (or strictly speaking double and add)
function *(n::U, P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem, U <: Integer}
    if isIdentity(P) || n == 1
        return P
    end
    E = P.curve
    if n == 0
        return identity(E)
    end
    if n < 0
        P = -P
        n = -n
    end
    bit = U(1) << (ndigits(n, base = 2) - 1)
    Q = P
    bit >>= 1
    while bit != 0
        Q = Q+Q
        if (bit & n) != 0
            Q += P
        end
        bit >>=1
    end
    return Q
end

function *(P::TwiHesECElem{T,S}, n::U) where {T <: CRing, S <: CRingElem, U <: Integer}
    return n*P
end

function ==(P::TwiHesECElem{T,S}, Q::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    P1 = normalized(P)
    Q1 = normalized(Q)
    return sameCurve(P,Q) && P1.X == Q1.X && P1.Y == Q1.Y && P1.Z == Q1.Z
end

############################
#
#   Elliptic Curve Stuff
#
############################

function jInvariant(E::TwiHesEC{T,S}) where {T <: CRing, S <: CRingElem}
    a = E.a
    d = E.d
    over = -(216*a*d+d^4)^3
    under = a*(27*a-d^3)^3
    return over*inv(under)
end

############################
#
#        Isogenies
#
############################

struct Isogeny
    domain::EllipticCurve
    codomain::EllipticCurve
    carryP::Function
    carryPfallback::Function
end

# Constructor
function Isogeny(E::EllipticCurve, P::EllipticCurveElem)
    if isIdentity(P)
        error("$P is the identity!")
    end
    if isIdentity(2*P)
        return deg2_isogeny(E,P)
    end
    if isIdentity(3*P)
        return deg3_isogeny(E,P)
    end
    error("$P is not of order 2 or 3")
end

#Later optimization: Normalize is only needed for computing the D of the isogenous
#curve. This is not necessary until the very last step. Maybe we can skip it?
#NB! Then formulas change, so be a bit careful (this one relies on Z = 1).
#Reasoning for the last part: Thm 1 from Dang and Moody says that order(P) = 3 <=> at one coord of P is zero
function deg2_isogeny(E::TwiHesEC{T,S}, P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    P1 = normalized(P)
    a = E.a
    d = E.d
    F = E.parent
    A = a^2
    s = P1.X
    t = P1.Y
    D = (-d+6*inv(s*t))*inv(s)
    E1 = TwiHesEC(F, A, D)
    phi(X, Y, Z) = (X*(X^2*t - s^2*Y*Z), Y*(Z^2*s*t - X*Y), Z*(Y^2*s - t^2*X*Z))
    phi2(X, Y, Z) =(X*(X*Z - Y^2*s*t), Y*(Y*Z*t^2 - X^2*a*s), Z*(X*Y*a*s^2 - Z^2*t))
    return Isogeny(E, E1, phi, phi2)
end

#Two different isogenies based on which subgroup we are looking at
function deg3_isogeny(E::TwiHesEC{T,S}, P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    if iszero(P.X)
        return deg3_isogeny_c1(E, P)
    else
        return deg3_isogeny_c2(E, P)
    end
end

#Case 1: {(0:-1:1),(0:-w:1),(0:-w^2:1)}, with w^3 = 1, w != 1
function deg3_isogeny_c1(E::TwiHesEC{T,S}, P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    a = E.a
    d = E.d
    F = E.parent
    A = d^3-27*a
    D = 3*d
    E1 = TwiHesEC(F, A, D)
    w = -normalized(P).Y
    phi(X, Y, Z) = (X*Y*Z, a*X^3 + w^2*Y^3 + w*Z^3, a*X^3 + w*Y^3 + w^2*Z^3)
    return Isogeny(E, E1, phi, phi)
end

#Case 2: {(0:-1:1),(1:-c:0),(1:0:-c)}, with c^3 = a
function deg3_isogeny_c2(E::TwiHesEC{T,S}, P::TwiHesECElem{T,S}) where {T <: CRing, S <: CRingElem}
    P1 = normalized(P)
    a = E.a
    d = E.d
    F = E.parent
    if iszero(P1.Z)
        c = -P1.Y
    else
        c = inv(-P1.X)
    end
    A = d^2*c + 3*d*c^2 + 9*a
    D = d + 6*c
    E1 = TwiHesEC(F, A, D)
    phi(X, Y, Z) = (X*Y*Z, c^2*X^2*Z + c*X*Y^2 + Y*Z^2, c^2*X^2*Y + c*X*Z^2 + Y^2*Z)
    return Isogeny(E, E1, phi, phi)
end

###### Printing ######
function show(io::IO, phi::Isogeny)
    print(io, "Isogeny φ : E_1 -> E_2 where \nE_1 = $(phi.domain)\nE_2 = $(phi.codomain)")
end

###### Go brrr ######
# Evaluate a point
function (phi::Isogeny)(P::EllipticCurveElem)
    if P.curve != phi.domain
        error("Point $P is not from the right curve")
    end
    (X, Y, Z) = phi.carryP(P.X,P.Y,P.Z)
    if !isZeroCoordinate(X, Y, Z)
        return phi.codomain(X, Y, Z)
    else
        (X, Y, Z) = phi.carryPfallback(P.X,P.Y,P.Z)
        return phi.codomain(X, Y, Z)
    end
end

############################
#
#           Other
#
############################
function allPointsNaive(E::TwiHesEC{T,S}) where {T <: CRing, S <: CRingElem}
    points = TwiHesECElem{T,S}[]
    Z = one(E.parent)
    for X in allElements(E.parent)
        for Y in allElements(E.parent)
            if validPoint(X, Y, Z, E)
                push!(points,E(X, Y, Z))
            end
        end
    end
    Z = zero(E.parent)
    X = one(E.parent)
    for Y in allElements(E.parent)
        if validPoint(X, Y, Z, E)
            push!(points,E(X, Y, Z))
        end
    end
    return points
end

end


text:
using OrdinaryDiffEqExtendedTests, OrdinaryDiffEq, Test

const CACHE_TEST_ALGS = [Euler(),Midpoint(),RK4(),SSPRK22(),SSPRK33(),SSPRK53(),
  SSPRK63(),SSPRK73(),SSPRK83(),SSPRK432(),SSPRK932(),SSPRK54(),SSPRK104(),CarpenterKennedy2N54(),
  BS3(),BS5(),DP5(),DP5Threaded(),DP8(),Feagin10(),Feagin12(),Feagin14(),TanYam7(),
  Tsit5(),TsitPap8(),Vern6(),Vern7(),Vern8(),Vern9(),OwrenZen3(),OwrenZen4(),OwrenZen5()]

#@time @testset "Feagin Tests" begin include("ode_feagin_tests.jl") end
@time @testset "Units Tests" begin include("units_tests.jl") end
@time @testset "Time Derivative Tests" begin include("time_derivative_test.jl") end
@time @testset "Jacobian Tests" begin include("jacobian_tests.jl") end
@time @testset "Number Type Tests" begin include("ode_numbertype_tests.jl") end
@time @testset "Ndim Complex Tests" begin include("ode_ndim_complex_tests.jl") end
@time @testset "In-Place Tests" begin include("ode_inplace_tests.jl") end
@time @testset "Event Detection Tests" begin include("event_detection_tests.jl") end
@time @testset "Autodiff Event Tests" begin include("autodiff_events.jl") end
@time @testset "Adaptive Tests" begin include("ode_adaptive_tests.jl") end
@time @testset "PSOS and Energy Conserve" begin include("psos_and_energy_conservation.jl") end
@time @testset "Init dt vs dorpri tests" begin include("init_dt_vs_dopri_tests.jl") end
@time @testset "Unrolled Tests" begin include("ode_unrolled_comparison_tests.jl") end
@time @testset "Alg Events Tests" begin include("alg_events_tests.jl") end


text:
__precompile__()
module JuliaDB

using IndexedTables, Dagger, NamedTuples, OnlineStats

import Base: collect, select, join
import IndexedTables: NextTable, table, NDSparse, ndsparse, Tup, groupjoin
import TextParse: csvread
import IndexedTables: Table
import Dagger: compute, distribute, free!, gather, load, save
using DataValues

# re-export
export IndexedTable, AbstractNDSparse, NDSparse, NextTable, Columns, colnames,
       table, ndsparse, compute, groupby, summarize, groupreduce, groupjoin,
       ColDict, insertafter!, insertbefore!, @cols, setcol, pushcol,
       popcol, insertcol, insertcolafter, insertcolbefore,
       renamecol, NA, dropna, flatten, ML, All, Not, Between, Keys

include("util.jl")
include("serialize.jl")
include("interval.jl")
include("table.jl")
include("ndsparse.jl")
include("reshape.jl")

# equality

function (==)(x::DDataset, y::Union{Dataset, DDataset})
    y1 = distribute(y, length.(domainchunks(rows(x))))
    res = delayed(==, get_result=true).(x.chunks, y1.chunks)
    all(collect(delayed((xs...) -> [xs...])(res...)))
end
function (==)(x::DDataset, y::Dataset)
    collect(x) == y
end
(==)(x::Dataset, y::DDataset) = y == x

include("iteration.jl")
include("sort.jl")

include("io.jl")
include("printing.jl")

include("indexing.jl")
include("selection.jl")
include("reduce.jl")
include("flatten.jl")
include("join.jl")

include("diagnostics.jl")
include("recipes.jl")
include("ml.jl")

end # module


text:
using Images
using VideoIO

f = VideoIO.openvideo("test.mp4")

# img = read(f,Image)


text:
function to_space(obj::Embedding, arg0::Point)
    return jcall(obj, "toSpace", Point, (Point,), arg0)
end

function to_sub_space(obj::Embedding, arg0::Point)
    return jcall(obj, "toSubSpace", Point, (Point,), arg0)
end


text:
"return `x+1` with probability `p` or `x-1` with `1-p`."
function update(x::Real,p::Real=0.5)
    ifelse(rand() < p, x+one(x), x-one(x))
end

max_stars_repo_name: isuke/JuliaBook-Samples

max_issues_repo_name: isuke/JuliaBook-Samples

max_forks_repo_name: isuke/JuliaBook-Samples


text:
###############
# H-OBCA: Hierarchical Optimization-based Collision Avoidance - a path planner for autonomous parking
# Copyright (C) 2018
# Alexander LINIGER [liniger@control.ee.ethz.ch; Automatic Control Lab, ETH Zurich]
# Xiaojing ZHANG [xiaojing.zhang@berkeley.edu; MPC Lab, UC Berkeley]
# Atsushi SAKAI [atsushisakai@global.komatsu; Komatsu Ltd / MPC Lab]
#
# 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 3 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/>.
###############
# The paper describing the theory can be found here:
# 	X. Zhang, A. Liniger and F. Borrelli; "Optimization-Based Collision Avoidance"; Technical Report, 2017, [https://arxiv.org/abs/1711.03449]
#   X. Zhang, A. Liniger, A. Sakai and F. Borrelli; "Autonomous  Parking  using  Optimization-Based  Collision  Avoidance"; Technical Report, 2018 [add URL]
###############

println("loading utility functions ....")
include("modules/setup.jl")
close("all")

# define scenarios
#scenario = "parallel"
scenario = "backwards"

# fixed or variable sampling time 1/0
fixTime = 0		# default: 0 (variable time steps)

#### problem parameters ####
TsPF = 0.05
if scenario == "backwards"
	# nominal sampling time
	sampleN = 3 		# down-sampling from Hybrid A* to OBCA
	if fixTime == 1 	#
		Ts = 0.65/3*sampleN		# 0.65/3 must be compatible with motion resolution of Hybrid A* algorithm
	else
		Ts = 0.6/3*sampleN 		# 0.6/3 must be compatible with motion resolution of Hybrid A* algorithm
	end
else
	sampleN = 3
	if fixTime == 1
		Ts = 0.95/3*sampleN		# 0.95/3 must be compatible with motion resolution of Hybrid A* algorithm
	else
		Ts = 0.9/3*sampleN		# 0.9/3 must be compatible with motion resolution of Hybrid A* algorithm
	end
end

# wheelbase
L  = 2.7

motionStep = 0.1	# step length of Hybrid A*",

# "nominal" shape of ego/controlled car,  ego object is later rotated around the car center
# center of rear wheel axis is reference point
# size of car is: (x_upper + x_lower) + (y_upper + y_lower)
#	   [x_upper, y_upper, -x_lower, -y_lower ]
ego  = [ 3.7   , 1      ,  1      ,  1       ]

##### define obstacles; for simplicity, only polyhedral obstacles are supported at this point
# obstacles are defined by vertices, which are assumed to be enumerated in clock-wise direction
# the first vertex must appear at the end of the list

# for plotting
nObPlot =  3 		# number of obstacles
vObPlot = [4 4 4]	# number of vertices of each obstacle, vector of dimenion nOb

# obstacle representation for optimization problem
nOb =  3 			# number of obstacles
vOb = [3 3 2]		# number of vertices of each obstacle, vector of dimenion nOb
vObMPC = vOb-1		# adjustment for optimizaton problem

if scenario == "backwards"
	println("Start Reverse Parking")
elseif scenario == "parallel"
	println("Start Parallel Parking")
else
	println("ERROR: please specify parking scenario")
end

# build environment
if scenario == "backwards"
	# obstacles for backwards
	#     	[ 	[[obst1_x1;obst1_y1],[obst1_x2;obst1_y2],[obst1_x3;obst1_y4],...,[obst1_x1;obst1_y1]]    , 		[[obst2_x1;obst2_y1],[obst2_x2;obst2_y2],[obst2_x3;obst2_y4],...,[obst2_x1;obst2_y1]]     ,     ...   ]
	lObPlot = [   [ [-20;5], [-1.3;5], [-1.3;-5], [-20;-5], [-20;5] ]  ,
	 	  [ [1.3;5], [20;5], [20;-5], [1.3;-5], [1.3;5] ] ,
		  [ [-20;15], [20;15], [20;11], [-20,11], [-20;15] ]		]		#vetices given in CLOCK-WISE direction

    # for optimization problem
	lOb = [   [ [-20;5], [-1.3;5], [-1.3;-5]]  ,
	 	  [ [1.3;-5] , [1.3;5] , [20;5] ] ,
		  [ [20;11], [-20;11]]	]		#vetices given in CLOCK-WISE direction

	# final state
	xF = [ 0 1.3 pi/2 0]

	# build obstacles for Hybrid A* algorithm
	ox = Float64[]
	oy = Float64[]
	# obstacle 1
	for i = -12:0.1:-1.3
	    push!(ox, Float64(i))
	    push!(oy, 5.0)
	end
	for i in -2:5
	    push!(ox, -1.3)
	    push!(oy, Float64(i))
	end
	# obstacle 2
	for i in -2:5
	    push!(ox, 1.3)
	    push!(oy, Float64(i))
	end
	for i = 1.3:0.1:12
	    push!(ox, Float64(i))
	    push!(oy, 5.0)
	end
	# obstacle 3
	for i = -12:12
	    push!(ox, Float64(i))
	    push!(oy, 11.0)
	end

elseif scenario == "parallel"
	# obstacles for backwards
	#     	[ 	[[obst1_x1;obst1_y1],[obst1_x2;obst1_y2],[obst1_x3;obst1_y4],...,[obst1_x1;obst1_y1]]    , 		[[obst2_x1;obst2_y1],[obst2_x2;obst2_y2],[obst2_x3;obst2_y4],...,[obst2_x1;obst2_y1]]     ,     ...   ]
	lObPlot = [   [ [-15;5], [-3;5], [-3;0], [-15;0], [-15;5] ]  ,
		 	  [ [3;5], [15;5], [15;0], [3;0], [3;5] ] ,
			  [ [-3;0], [-3;2.5], [3;2.5], [3,0], [-3;0] ]		]		#vetices given in CLOCK-WISE direction
	# obstacle representation for optimization problem
  	lOb = [   [ [-20;5], [-3.;5], [-3.;0]]  ,
	 	  [ [3.;0] , [3.;5] , [20;5] ] ,
		  [ [-3;2.5], [ 3;2.5]]]		#vetices given in CLOCK-WISE direction

	# final state
	xF = [-L/2 4 0 0]

	# obstacles for Hybrid A* algorithms
	ox = Float64[]
	oy = Float64[]
	# obstacle 1
	for i in -12:0.1: -3.
		push!(ox,Float64(i))
		push!(oy,5.0)
	end
	for i in -2  : 5
		push!(ox,-3.0)
		push!(oy,Float64(i))
	end
	# obstacle 2
	for i in -3 : 3
		push!(ox,Float64(i))
		push!(oy,2.5)
	end
	# obstacle 3
	for i in -2 : 5
		push!(ox,3.0)
		push!(oy,Float64(i))
	end
	for i in 3 :0.1: 12
		push!(ox,Float64(i))
		push!(oy,5.0)
	end
	# obstacle 4
	for i in -12 : 12
		push!(ox,Float64(i))
		push!(oy,11.0)
	end
end

#           [x_lower, x_upper, -y_lower,   y_upper  ]
XYbounds =  [ -15   , 15      ,  1      ,  10       ]

# set initial state
x0 = [-6  9.5   0.0    0.0]
# x0 = [9  7   0.0    0.0]

# call Hybrid A*
tic()
rx, ry, ryaw = hybrid_a_star.calc_hybrid_astar_path(x0[1], x0[2], x0[3], xF[1], xF[2], xF[3], ox, oy, hybrid_a_star.XY_GRID_RESOLUTION, hybrid_a_star.YAW_GRID_RESOLUTION, hybrid_a_star.OB_MAP_RESOLUTION)
timeHybAstar = toq();

### extract (smooth) velocity profile from Hybrid A* solution ####
rv = zeros(length(rx),1)
for i=1:length(rx)
	if i < length(rx)
		rv[i] = (rx[i+1] - rx[i])/(Ts/sampleN)*cos(ryaw[i]) + (ry[i+1]-ry[i])/(Ts/sampleN)*sin(ryaw[i])
	else
		rv[i] = 0
	end
end
### Smoothen velocity 0.3 m/s^2 max acceleration ###
v,a = veloSmooth(rv,0.3,Ts/sampleN)
### compute steering angle ###
delta = atan.(diff(ryaw)*L/motionStep.*sign.(v[1:end-1]));

### Down-sample for Warmstart ##########
rx_sampled = rx[1:sampleN:end]
ry_sampled = ry[1:sampleN:end]
ryaw_sampled = ryaw[1:sampleN:end]
v_sampled = v[1:sampleN:end]

a_sampled = a[1:sampleN:end]
delta_sampled = delta[1:sampleN:end]

## initialize warm start solution
xWS = [rx_sampled ry_sampled ryaw_sampled v_sampled]
uWS = [delta_sampled a_sampled]

### solve OBCA step ###
N = length(rx_sampled)-1
AOb, bOb = obstHrep(nOb, vOb, lOb) 	# obtain H-representation of obstacles
xp10, up10, scaleTime10, exitflag10, time10, lp10, np10 = ParkingSignedDist(x0,xF,N,Ts,L,ego,XYbounds,nOb,vObMPC,AOb,bOb,fixTime,xWS,uWS)

### plot H-OBCA solution ###
if exitflag10==1
	println("H-OBCA successfully completed.")
	figure(1)
	hold(1)
	plot(xp10[1,:],xp10[2,:],"b")

	plotTraj(xp10',up10',length(rx_sampled)-1,ego,L,nObPlot,vObPlot,lObPlot,"Trajectory generated by H-OBCA",1)
else
	println("  WARNING: Problem could not be solved.")
end

### comparison with Hybrid A* ###
figure(2)
title("Trajectory Comparison")
hold(1)
for j = 1 : nObPlot # plot obstacles
	for k = 1 : vObPlot[j]
		plot([lObPlot[j][k][1],lObPlot[j][k+1][1]] , [lObPlot[j][k][2],lObPlot[j][k+1][2]] ,"k")
	end
end
plot(xp10[1,:],xp10[2,:], "-b",  label="H-OBCA")
plot(rx, ry, "--r",  label="Hybrid A*")
plot(x0[1],x0[2],"ob")
plot(xF[1],xF[2],"ob")
legend()
axis("equal")

totTime = timeHybAstar+time10	# total execution time of H-OBCA
println("Total run time: " , totTime, " s")
println("  Hybrid A* time: ", timeHybAstar, " s")
println("  optimization (OBCA) time: ", time10, " s")

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
module CommunityServiceUNAM

# package code goes here

end # module

max_stars_repo_name: miguelraz/CommunityServiceUNAM

max_issues_repo_name: miguelraz/CommunityServiceUNAM

max_forks_repo_name: miguelraz/CommunityServiceUNAM


text:
include("RedBlack.jl")
include("GeometricalWeights.jl")
include("Painting.jl")
include("PotentialCalculationSetup/PotentialCalculationSetup.jl")

include("SuccessiveOverRelaxation/SOR.jl")
include("Refinement.jl")


text:
# ==========================================================================
# Unprocessed Solution for the dual problem

immutable DualSolution
  minimizer::Vector{Float64}
  bias::Float64
  minimum::Float64
  iterations::Int
  isconverged::Bool
end

function Base.show(io::IO, r::DualSolution)
  println(io, typeof(r))
  _printvariable(io, 14, "minimizer()", minimizer(r))
  _printvariable(io, 14, "bias()", bias(r))
  _printvariable(io, 14, "minimum()", minimum(r))
  _printvariable(io, 14, "iterations()", iterations(r))
  _printvariable(io, 14, "isconverged()", isconverged(r), newline = false)
end

coef(s::DualSolution) = s.minimizer
minimizer(s::DualSolution) = s.minimizer
Base.minimum(s::DualSolution) = s.minimum
bias(s::DualSolution) = s.bias
iterations(s::DualSolution) = s.iterations
isconverged(s::DualSolution) = s.isconverged


text:
module problem114

description = """
Counting block combinations I
A row measuring seven units in length has red blocks with a minimum length of three units placed on it,
 such that any two red blocks (which are allowed to be different lengths) are separated by at least one black square.
There are exactly seventeen ways of doing this.
How many ways can a row measuring fifty units in length be filled?
"""

function block_combinations(n :: Integer, memo = zeros(Int, n))
  # trivial cases
  if n < 1 return 1 end

  # Count the total with each block size, and block position,
  # recursively counting the ways to fill the remaining space at the end of the row.
  if memo[n] == 0
    total = 1
    for blocksize = 3:n, blockpos = 1:(n-blocksize+1)
      total += block_combinations(n - blockpos - blocksize, memo)
    end
    memo[n] = total
  end
  return memo[n]
end

using Base.Test

@test block_combinations(3) == 2 #empty or 1 block
@test block_combinations(4) == 4 #empty or length 3 at pos 1, or length 3 at pos 2, or length 4 at pos 1
@test block_combinations(7) == 17

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
using LinearAlgebra

mutable struct Node
    index::Int64
    centroid::Vector{Float64}
    iFaces::Vector{Int64}
    iCells::Vector{Int64}
    Node() = new(0, [], [], [])
end

mutable struct Face
    index::Int64
    iOwner::Int64
    iNeighbour::Int64
    iNodes::Vector{Int64}
    centroid::Vector{Float64}
    area::Float64
    Sf::Vector{Float64}
    n::Vector{Float64}
    dCF::Vector{Float64}
    eCF::Vector{Float64}
    dCf::Vector{Float64}
    Tf::Vector{Float64}
    geoDiff::Float64
    wallDist::Float64
    gf::Float64
    iOwnerNeighbourCoeff::Int64
    iNeighbourOwnerCoeff::Int64
    Face() = new(0, 0, 0, [], [], 0.0, [], [], [], [], [], [], 0.0, 0.0, 0.0, 0, 0)
end

mutable struct Cell
    index::Int64
    iNeighbours::Vector{Int64}
    iNodes::Vector{Int64}
    iFaces::Vector{Int64}
    volume::Float64
    centroid::Vector{Float64}
    faceSigns::Vector{Int64}
    nNeighbours::Int64
    Cell() = new(0, [], [], [], 0.0, [], [], 0)
end

mutable struct Boundary
    index::Int64
    name::String
    type::String
    startFace::Int64
    nFaces::Int64
    bcType::String
    Boundary() = new(0, "", "", 0, 0, "")
end

mutable struct Mesh
    nodes::Vector{Node}
    faces::Vector{Face}
    cells::Vector{Cell}
    boundaries::Vector{Boundary}
    nNodes::Int64
    nFaces::Int64
    nInternalFaces::Int64
    nBoundaryFaces::Int64
    nCells::Int64
    nBoundaries::Int64
end

function nodeIndexToCoord(iNodes::Vector{Int64}, nodes::Vector{Vector{Float64}})
    nodeCoords = Vector{Float64}[]
    for iNode in iNodes
        push!(nodeCoords, nodes[iNode])
    end
    return nodeCoords
end

function faceIndexToSf(iFaces::Vector{Int64}, faceVec::Vector{Face})
    facesSf = Vector{Float64}[]
    for iFace in iFaces
        push!(facesSf, faceVec[iFace].Sf)
    end
    return facesSf
end

function faceIndexToCentroid(iFaces::Vector{Int64}, faceVec::Vector{Face})
    facesCentroid = Vector{Float64}[]
    for iFace in iFaces
        push!(facesCentroid, faceVec[iFace].centroid)
    end
    return facesCentroid
end

function geometricCenter(points::Vector{Vector{Float64}})
    nPoints = size(points, 1)
    gc = [0, 0, 0]
    for point in points
        gc = gc .+ point
    end
    return gc / nPoints
end

function triangleSurfaceVec(points::Vector{Vector{Float64}})
    side1 = points[2] .- points[1]
    side2 = points[3] .- points[1]
    return cross(side1, side2) / 2
end

function polygonCentroidSurfaceVec(points::Vector{Vector{Float64}})
    nPoints = size(points, 1)
    gc = geometricCenter(points)

    centroid = [0.0, 0.0, 0.0]
    surfaceVec = [0.0, 0.0, 0.0]

    if nPoints < 3
        throw(ErrorException("Error calculating surface centroid: The number of points is less than 3"))
    elseif nPoints == 3
        centroid = gc
        surfaceVec = triangleSurfaceVec(points)
    else
        for iPoint in 1:nPoints
            if iPoint < nPoints
                subTriPoints = [gc, points[iPoint], points[iPoint + 1]]
            else
                subTriPoints = [gc, points[iPoint], points[1]]
            end

            subTrisurfaceVec = triangleSurfaceVec(subTriPoints)
            subTriCentroid = geometricCenter(subTriPoints)

            centroid = centroid .+ (subTriCentroid * norm(subTrisurfaceVec))
            surfaceVec = surfaceVec .+ subTrisurfaceVec

        end
        centroid /= norm(surfaceVec)
    end
    return centroid, surfaceVec
end

function polyhedralCentroidVol(points, surfaceVecs, surfaceCentroids)
    gc = geometricCenter(points)
    nFaces = size(surfaceVecs, 1)

    vol = 0.0
    centroid = [0.0, 0.0, 0.0]

    for i in 1:nFaces
        geometriCenterToFaceVec = gc .- surfaceCentroids[i]
        subPyrVol = abs(dot(geometriCenterToFaceVec, surfaceVecs[i]) / 3)
        centroidSubPyr = 0.75 * surfaceCentroids[i] + 0.25 * gc

        vol += subPyrVol
        centroid = centroid .+ (centroidSubPyr * subPyrVol)
    end
    centroid /= vol
    return centroid, vol
end

function processOpenFoamMesh(points, faces, owner, neighbour, boundary)
    nNodes = length(points)
    nFaces = length(faces)
    nInternalFaces = length(neighbour)
    nBoundaryFaces = nFaces - nInternalFaces
    nBoundaries = length(boundary)
    maximum(owner) > maximum(neighbour) ? nCells = maximum(owner) : nCells = maximum(neighbour)

    nodeVec = [Node() for _ in 1:nNodes]
    faceVec = [Face() for _ in 1:nFaces]
    cellVec = [Cell() for _ in 1:nCells]
    boundaryVec = [Boundary() for _ in 1:nBoundaries]

    for iFace in 1:nFaces
        faceVec[iFace].index = iFace
        faceVec[iFace].iNodes = faces[iFace]
        faceVec[iFace].iOwner = owner[iFace]

        faceNodesCoord = nodeIndexToCoord(faces[iFace], points)
        faceVec[iFace].centroid, faceVec[iFace].Sf = polygonCentroidSurfaceVec(faceNodesCoord)
        faceVec[iFace].area = norm(faceVec[iFace].Sf)

        cellVec[owner[iFace]].index = owner[iFace]
        push!(cellVec[owner[iFace]].iFaces, iFace)
        push!(cellVec[owner[iFace]].faceSigns, 1)

        for iNode in faceVec[iFace].iNodes
            nodeVec[iNode].index = iNode
            nodeVec[iNode].centroid = points[iNode]
            push!(nodeVec[iNode].iFaces, iFace)
            !(owner[iFace] in nodeVec[iNode].iCells) && push!(nodeVec[iNode].iCells, owner[iFace])
            !(iNode in cellVec[owner[iFace]].iNodes) && push!(cellVec[owner[iFace]].iNodes, iNode)

            if iFace <= nInternalFaces
                !(neighbour[iFace] in nodeVec[iNode].iCells) && push!(nodeVec[iNode].iCells, neighbour[iFace])
                !(iNode in cellVec[neighbour[iFace]].iNodes) && push!(cellVec[neighbour[iFace]].iNodes, iNode)
            end
        end

        if iFace <= nInternalFaces
            faceVec[iFace].iNeighbour = neighbour[iFace]

            push!(cellVec[neighbour[iFace]].iFaces, iFace)
            push!(cellVec[neighbour[iFace]].faceSigns, -1)

            push!(cellVec[owner[iFace]].iNeighbours, neighbour[iFace])
            cellVec[owner[iFace]].nNeighbours += 1

            push!(cellVec[neighbour[iFace]].iNeighbours, owner[iFace])
            cellVec[neighbour[iFace]].nNeighbours += 1
        else
            faceVec[iFace].iNeighbour = -1
        end
    end

    for iCell in 1:nCells
        cellNodesCoord = nodeIndexToCoord(cellVec[iCell].iNodes, points)
        cellFacesSf = faceIndexToSf(cellVec[iCell].iFaces, faceVec)
        cellFacesCentroid = faceIndexToCentroid(cellVec[iCell].iFaces, faceVec)
        cellVec[iCell].centroid, cellVec[iCell].volume = polyhedralCentroidVol(cellNodesCoord, cellFacesSf, cellFacesCentroid)
    end

    for iBoundary in 1:nBoundaries
        boundaryVec[iBoundary].index = iBoundary
        boundaryVec[iBoundary].name = boundary[iBoundary]["name"]
        boundaryVec[iBoundary].type = boundary[iBoundary]["type"]
        boundaryVec[iBoundary].startFace = boundary[iBoundary]["startFace"]
        boundaryVec[iBoundary].nFaces = boundary[iBoundary]["nFaces"]
    end

    for face in faceVec
        face.dCf = face.centroid - cellVec[face.iOwner].centroid
        face.n = face.Sf/(norm(face.Sf))

        if face.iNeighbour != -1
            face.dCF = cellVec[face.iNeighbour].centroid - cellVec[face.iOwner].centroid
            dfn = cellVec[face.iNeighbour].centroid - face.centroid
            face.gf = dot(face.dCf, face.n) / (dot(face.dCf, face.n) + dot(dfn, face.n))
            face.wallDist = 0
        else
            face.dCF = face.dCf
            face.gf = 1
            face.wallDist = dot(face.dCf, face.n)
        end
        face.eCF = face.dCF / norm(face.dCF)
        Ef = face.area * face.eCF
        face.geoDiff = face.iNeighbour != -1 ? norm(Ef) / norm(face.dCF) : norm(Ef) / face.wallDist
        face.Tf = face.Sf - Ef
    end

    return Mesh(nodeVec, faceVec, cellVec, boundaryVec, nNodes,
        nFaces, nInternalFaces, nBoundaryFaces, nCells, nBoundaries)
end

function printNode(node)
    println("")
    println("Node $(node.index) = ")
    println("")
    println("                   index: ", node.index)
    println("                centroid: ", node.centroid)
    println("                  iFaces: ", node.iFaces)
    println("                  iCells: ", node.iCells)
end

function printFace(face)
    println("")
    println("Face $(face.index) = ")
    println("")
    println("                   index: ", face.index)
    println("                  iNodes: ", face.iNodes)
    println("                  iOwner: ", face.iOwner)
    println("              iNeighbour: ", face.iNeighbour)
    println("                centroid: ", face.centroid)
    println("                    area: ", face.area)
    println("                      Sf: ", face.Sf)
    println("                      CN: ", face.dCf)
    println("                 geodiff: ", face.geoDiff)
    println("                       T: ", face.Tf)
    println("                      gf: ", face.gf)
    println("                walldist: ", face.wallDist)
    println("    iOwnerNeighbourCoeff: ", face.iOwnerNeighbourCoeff)
    println("    iNeighbourOwnerCoeff: ", face.iNeighbourOwnerCoeff)
end

function printBoundary(boundary)
    println("")
    println("Boundary $(boundary.index) = ")
    println("")
    println("                   index: ", boundary.index)
    println("                    name: ", boundary.name)
    println("                    type: ", boundary.type)
    println("               startFace: ", boundary.startFace)
    println("          numberOfBFaces: ", boundary.nFaces)
end

function printCell(cell)
    println("")
    println("Cell $(cell.index) = ")
    println("")
    println("                   index: ", cell.index)
    println("                  iFaces: ", cell.iFaces)
    println("                  iNodes: ", cell.iNodes)
    println("             iNeighbours: ", cell.iNeighbours)
    println("      numberOfNeighbours: ", cell.nNeighbours)
    println("               faceSigns: ", cell.faceSigns)
    println("                  volume: ", cell.volume)
    println("                centroid: ", cell.centroid)
end

function printMesh(mesh)
    println("")
    println("Mesh = ")
    println("")
    println("                   nodes: ", typeof(mesh.nodes))
    println("           numberOfNodes: ", mesh.nNodes)
    println("                   faces: ", typeof(mesh.faces))
    println("           numberOfFaces: ", mesh.nFaces)
    println("                   cells: ", typeof(mesh.cells))
    println("           numberOfCells: ", mesh.nCells)
    println("              boundaries: ", typeof(mesh.boundaries))
    println("      numberOfBoundaries: ", mesh.nBoundaries)
    println("   numberOfInteriorFaces: ", mesh.nInternalFaces)
    println("   numberOfBoundaryFaces: ", mesh.nBoundaryFaces)
end


text:
using PSCOPF

function get_test_case()
    network = PSCOPF.Networks.Network()
    # Buses : BUS
    PSCOPF.Networks.add_new_bus!(network, "bus_1")
    PSCOPF.Networks.add_new_bus!(network, "bus_2")
    # Branches : BRANCH LIMIT
    PSCOPF.Networks.add_new_branch!(network, "branch_1_2", 500.);
    # PTDF : BRANCH BUS VALUE
    PSCOPF.Networks.add_ptdf_elt(network, "branch_1_2", "bus_1", 0.5)
    PSCOPF.Networks.add_ptdf_elt(network, "branch_1_2", "bus_2", -0.5)
    #Alternatively,
    # Generators
    #Limitables
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_0", PSCOPF.Networks.LIMITABLE,
                                            0., 0., #pmin, pmax : Not concerned ? min is always 0, max is the limitation
                                            0., 10., #start_cost, prop_cost : start cost is always 0 ?
                                            Dates.Second(0), Dates.Second(0)) #dmo, dp : always 0. ?
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "wind_2_0", PSCOPF.Networks.LIMITABLE,
                                            0., 0.,
                                            0., 11.,
                                            Dates.Second(0), Dates.Second(0))
    #Imposables
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "ccg_1_0", PSCOPF.Networks.IMPOSABLE,
                                            10., 200., #pmin, pmax
                                            45000., 30., #start_cost, prop_cost
                                            Dates.Second(4*3600), Dates.Second(15*60)) #dmo, dp
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "tac_2_0", PSCOPF.Networks.IMPOSABLE,
                                            10., 200.,
                                            12000., 120.,
                                            Dates.Second(30*60), Dates.Second(15*60))

    #initial generators state
    generators_init_state = SortedDict(
                    "ccg_1_0" => PSCOPF.ON,
                    "tac_2_0" => PSCOPF.OFF,
                )

    # Uncertainties : ECH NAME TS S VALUE
    uncertainties = PSCOPF.Uncertainties()
    #Alternatively, PSCOPF.add_uncertainty!(uncertainties, ech, nodal_injection_name, ts, scenario_name, value)

    # Timesteps
    TS = [DateTime("2015-01-01T11:00:00"),
            DateTime("2015-01-01T11:15:00"),
            DateTime("2015-01-01T11:30:00"),
            DateTime("2015-01-01T11:45:00")]

    mode = PSCOPF.PSCOPF_MODE_1
    ECH = PSCOPF.generate_ech(network, TS, mode)

    sequence = PSCOPF.generate_sequence(network, TS, ECH, mode)

    exec_context = PSCOPF.PSCOPFContext(network, TS, ECH, mode, generators_init_state, uncertainties, nothing)
    PSCOPF.add_schedule!(exec_context, PSCOPF.Schedule(PSCOPF.Market(), ECH[1]))
    PSCOPF.add_schedule!(exec_context, PSCOPF.Schedule(PSCOPF.TSO(), ECH[1]))
    PSCOPF.run!(exec_context, sequence)

    market_schedule = PSCOPF.safeget_last_market_schedule(exec_context)
    tso_schedule = PSCOPF.safeget_last_tso_schedule(exec_context)

end
#=
    Description of all optimization problems
=#

# ech : value
# TS : list
# SCENARIO : list
# commitment :
function energy_market(ech, TS, SCENARIO, commitment)
#=
    availability is given by DMO/DP and ECH

    minimize production cost + fixed cost
    s.c.
        EOD[ts, s] for all ts in TS, s in SCENARIO
    output :
        - commitment market for DMO=ECH or commitment[s] : this is because market realization cannot be forecasted
        - P for DP=ECH or P[s]
=#
end

# ech : value
# TS : list
# SCENARIO : list
# commitment :
function tso_out_operationnal_window(ech, TS, SCENARIO, commitment)
    #=
        imposition are related to DP/DMO, cost are computed by delta = |P[s]-P_imp|

        1- sum delta_up+delta_down
        2- minimize cost to have

        minimize 1- or 2- with 1- in constraints
        s.c.
            energy_market result EOD OK
            line constraints
            equilibrated imposition (computed how ? ) sum delta = 0

        output :
            - commitment tso for DMO=ECH or commitment[s]
            - P for DP=ECH, P[s] tso
    =#
end

# ech : value
# TS : list
# SCENARIO : value or list
# commitment :
# question : one scenario to be sure it is feasible (average, other ?), or multi-scenario with delta allowed
function entering_operationnal_window(ech, TS, SCENARIO, commitment)
    #=
    availability is given by DMO/DP and ECH

    minimize production cost
    s.c.
        EOD[ts] for all ts in TS

    output :
        - commitment, same value for all s in SCENARIO, market
        - P, same value for all s in SCENARIO, market
    =#
end

# ech : value
# TS : list
# SCENARIO : list
# commitment :
function tso_in_operationnal_window(ech, TS, SCENARIO, commitment)
    #=

    =#
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    int[] foo = new int[1];

    void m() {
	final Condition c = new Condition();
	int y = ({L /c T}) 42;
	foo[0] = y;
    }

    public static void main(String[] args) {
	new C().m();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
module RulesBernoulliPTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "rules:Bernoulli:p" begin

    @testset "Belief Propagation: (m_out::PointMass)" begin

        @test_rules [ with_float_conversions = true ] Bernoulli(:p, Marginalisation) [
            (input = (m_out = PointMass(1.0),), output = Beta(2.0, 1.0)),
            (input = (m_out = PointMass(0.2),), output = Beta(12/10, 9/5)),
        ]

    end

    @testset "Variational Message Passing: (q_out::PointMass)" begin

        @test_rules [ with_float_conversions = true ] Bernoulli(:p, Marginalisation) [
            (input = (q_out = Bernoulli(1.0),), output = Beta(2.0, 1.0)),
            (input = (q_out = Bernoulli(0.3),), output = Beta(13/10, 17/10)),
        ]

    end
end
end


text:
mutable struct MMVMNFAModel
    name::String
    pivec::Matrix{Float64}
    mu::Matrix{Float64}
    beta::Matrix{Float64} #Stochastic mean component
    B::Union{Array{Float64},Array{ComplexF64}}
    D::Array{Float64}
    Ψ::Union{Matrix{Float64},Nothing} #Parameters of the scaling density function
    g::Int64
    q::Int64
    tau::Matrix{Float64}
    w_esteps::Union{Array{Float64},Nothing} #Esteps for the scaling density function
    logL::Float64
    BIC::Float64
    status::String
    niter::Int64
    logL_history::Vector{Float64}
end

mutable struct BUpdateQuantities
    Lambda_cs::Matrix{Float64}
    n_param::Vector{Int64}
    q_scores::Matrix{Float64}
    q_up::Int64
    q_change::Bool
    B_update::Array{ComplexF64}
end

mutable struct MtFAQuantities
    fnu::Any
    optim_min::Any #Define the optim min class here... and union with nothing
end

mutable struct MSLFAQuantities
    integrand::Any
    int_res::Float64
    c::Vector{Float64} #possible vector??
    integrand_estep::Any
    int_res_estep::Vector{Float64}
end

mutable struct MCNFAQuantities
    d1::Vector{Float64}
    d2::Vector{Float64}
    kappa_ij::Vector{Float64}
end

mutable struct MGHFAQuantities
    zeta_bar::Float64
    xi_bar::Float64
    rho_bar::Float64
    flam::Any
    fomega::Any
    min_lam::Any #Define the optim min class here... and union with nothing
    min_omega::Any #Define the optim min class here... and union with nothing
    delta_b::Float64
    omega::Float64
    lambda::Float64
    t1::Vector{Float64}
    t2::Vector{Float64}
    integrand::Any
    t3::Vector{Float64}
end

mutable struct MBSFAQuantities
    xi_bar::Float64
    rho_bar::Float64
    a_is::Float64
    delta_b::Float64
    Theta::Vector{Float64}
    Lambda::Float64
    a::Float64
    b1::Vector{Float64}
    b2::Vector{Float64}
end

mutable struct MLFAQuantities
    xi_bar::Float64
    a::Float64
    delta_b::Float64
    numerator::Vector{Float64}
    py::Vector{Float64}
    rat::Vector{Float64}
    b1::Vector{Float64}
    b2::Vector{Float64}
    b3::Vector{Float64}
end

mutable struct MMVMNFAMutableQuantities
    BIC_best::Float64
    best_model::Union{Nothing, MMVMNFAModel}
    g::Int64
    q::Int64
    k_starts::Array{Float64,2}
    r_starts::Array{Float64,2}
    starts::Array{Float64,2}
    ind::Int64
    init_para_error::Union{String,Bool}
    indices::Any
    Si::Matrix{Float64}
    Di_sqrt::Matrix{Float64}
    inv_Di_sqrt::Matrix{Float64}
    lambda::Any
    eig_order::Any
    lambda_vals::Any
    sigma2::Float64
    maxinit::Int64
    sigma::Matrix{Float64}
    delta::Vector{Float64}
    Fji::Matrix{Float64}
    Fjmax::Matrix{Float64}
    n_i::Matrix{Float64}
    Vtilde::Array{Float64}
    Lambda::Matrix{Float64}
    Ymu::Array{Float64, 2}
    Vi::Array{Float64, 2}
    D_is::Any
    V_eig_temp::Any
    V_eig::Any
    qi_break::Bool
    qi::Int64
    Uqi::Any
    Ri::Any
    C::Any
    Us::Any
    V::Any
    bvec::Any
    psitilde::Array{Float64,2}
    c::Any
    v::Any
    b::Any
    a::Any
    omega::Any
    ratio::Any
    eigenvecs::Array{Any,1}
    logL_old::Float64
    d_model::Int64
    tau_rho_xi_ij::Matrix{Float64}
    tau_rho_ij::Matrix{Float64}
    tr_ij::Matrix{Float64}
    tx_ij::Matrix{Float64}
    xi_bar_i::Float64
    rho_bar_i::Float64
    beta_matrix::Matrix{Float64}
    mq::Union{Nothing,MtFAQuantities,MSLFAQuantities,MCNFAQuantities,MGHFAQuantities,MBSFAQuantities,MLFAQuantities}
    tau_xi_ij::Matrix{Float64}
    qup::Union{Nothing,BUpdateQuantities}
end

struct MMVMNFAStaticQuantities
    Y::Array{Float64, 2}
    gmin::Int64
    gmax::Int64
    qmin::Int64
    qmax::Int64
    eta::Float64
    itmax::Int64
    nkmeans::Int64
    nrandom::Int64
    psi_initial::Union{Nothing, Vector{Float64}}
    update_q::Bool
    tol::Float64
    conv_measure::String
    n::Int64
    p::Int64
    n_it::Int64
    n_comb::Int64
    combs::Union{Nothing, Matrix{Tuple{Int64, Int64}}}
    ecm_fit!::Any
    logfY!::Any
    esteps!::Any
    msteps::Any
    n_esteps::Union{Int64, Nothing}
    logL_calc!::Any
    Mstep_calc!::Any
end

max_stars_repo_name: john-c-davey/FactorMixtures

max_issues_repo_name: john-c-davey/FactorMixtures

max_forks_repo_name: john-c-davey/FactorMixtures


text:
function det_by_lu(y, x, N)
    y[1] = 1.

    for k = 1:N
        y[1] *= x[k,k]
        for i = k+1:N
            x[i,k] /= x[k,k]
	    end
        for j = k+1:N
            for i = k+1:N
                x[i,j] -= x[i,k] * x[k,j]
            end
        end
    end
end

nb_loops = []

function run_julia(y,A,B,N)
    loops = max(10000000 // (N*N), 1)
    append!(nb_loops, loops)
    print(loops, '\n')
    for l in 1:loops
        B[:,:] = A
        det_by_lu(y, B, N)
    end
end

y = [0.0]

times = []

N=5
A = rand(N,N)
B = zeros(N,N)
# Warm up
@time run_julia(y,A,B,N)

N=5
A = rand(N,N)
B = zeros(N,N)
append!(times, @elapsed run_julia(y,A,B,N))
N=10
A = rand(N,N)
B = zeros(N,N)
append!(times, @elapsed run_julia(y,A,B,N))
N=30
A = rand(N,N)
B = zeros(N,N)
append!(times, @elapsed run_julia(y,A,B,N))
N=100
A = rand(N,N)
B = zeros(N,N)
append!(times, @elapsed run_julia(y,A,B,N))
N=200
A = rand(N,N)
B = zeros(N,N)
append!(times, @elapsed run_julia(y,A,B,N))
N=300
A = rand(N,N)
B = zeros(N,N)
append!(times, @elapsed run_julia(y,A,B,N))
N=400
A = rand(N,N)
B = zeros(N,N)
append!(times, @elapsed run_julia(y,A,B,N))
N=600
A = rand(N,N)
B = zeros(N,N)
append!(times, @elapsed run_julia(y,A,B,N))
N=1000
A = rand(N,N)
B = zeros(N,N)
append!(times, @elapsed run_julia(y,A,B,N))

print("\ntime per loop=\n", times ./ nb_loops[2:end])
# [5, 10, 30, 100, 200, 300, 400, 600, 1000])

max_stars_repo_name: exabl/sandbox

max_issues_repo_name: exabl/sandbox

max_forks_repo_name: exabl/sandbox


text: add_gigasecond(date::DateTime) = date + Millisecond(10^12)

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
include("../src/RandomField.jl")
using .RandomField

# chose input parameters:
n 		= 1000		# number of grid points
n_real	= 10000		# number of realizations
l 		= 30.0		# real space box is [-l,l]
f 		= gaus_rf	# spectral density

var_y(ν) = 2*sqrt(2pi)*(l/pi)*f(ν)

# create structures (allocate memory)
Y = Vector{ComplexF64}(undef,n)	# fourier realization
Z = Vector{ComplexF64}(undef,n)	# real realization

M_est = zeros(ComplexF64,n) # est. mean
C_est = zeros(ComplexF64,n)	# est. stationary covariance
F_est = zeros(ComplexF64,n)	# est. spectral density
C_tru = Vector{ComplexF64}(undef,n) # true stationary covariance
F_tru = Vector{ComplexF64}(undef,n)	# true spectral density

T, V = create_grids(n,l)

# draw realization
draw_realization!(Y,var_y,V)
ift!(Z,Y,n,l)

# zeroth order moment
for i=1:n_real
	draw_realization!(Y,var_y,V)
	ift!(Z,Y,n,l)
	M_est .+= Z
end
M_est ./= n_real

# first order moments
for i=1:n_real
	draw_realization!(Y,var_y,V)
	ift!(Z,Y,n,l)
	@. F_est += sqrt(1/8pi)*(pi/l) * Y * conj(Y)
	@. C_est += 0.5 * Z * conj(Z[div(n,2)])
end
F_est ./= n_real
C_est ./= n_real

@. F_tru = f(V)
ift!(C_tru,F_tru,n,l)

# print to file:
out = open("../data/output.csv", "w")
print(out, "T, V, ReY, ReZ, ImZ, M_est, C_est, C_tru, F_est, F_tru\n")
for i in 1:n
	print(out,
		"$(T[i]),"*
		"$(V[i]),"*
		"$(real(Y[i])),"*
		"$(real(Z[i])),"*
		"$(imag(Z[i])),"*
		"$(real(M_est[i])),"*
		"$(real(C_est[i])),"*
		"$(real(C_tru[i])),"*
		"$(real(F_est[i])),"*
		"$(real(F_tru[i]))\n")
end
close(out)


text:
# ---
# title: 1014. Best Sightseeing Pair
# id: problem1014
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/best-sightseeing-pair/description/>
# hidden: true
# ---
#
# Given an array `A` of positive integers, `A[i]` represents the value of the
# `i`-th sightseeing spot, and two sightseeing spots `i` and `j` have distance
# `j - i` between them.
#
# The _score_  of a pair (`i < j`) of sightseeing spots is (`A[i] + A[j] + i -
# j)` : the sum of the values of the sightseeing spots, **minus** the distance
# between them.
#
# Return the maximum score of a pair of sightseeing spots.
#
#
#
# **Example 1:**
#
#
#
#     Input: [8,1,5,2,6]
#     Output: 11
#     Explanation: i = 0, j = 2, A[i] + A[j] + i - j = 8 + 5 + 0 - 2 = 11
#
#
#
#
# **Note:**
#
#   1. `2 <= A.length <= 50000`
#   2. `1 <= A[i] <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    AbstractOptimizationMethod

Abstract type used to implement an interface for generic optimization methods.
"""
abstract type AbstractOptimizationMethod{T} end


text:
using MPI
using P4est_wrapper
using Test

# Initialize MPI if not initialized yet
if !MPI.Initialized()
    MPI.Init()
end

#############################################################################
# Definition of data structures and function callbacks
#############################################################################

mpicomm = P4est_wrapper.P4EST_ENABLE_MPI ? MPI.COMM_WORLD : Cint(0)

# Dummy callback
dummy_callback( ::Ptr{p4est_t}, which_tree::p4est_topidx_t, quadrant::Ptr{p4est_quadrant_t}) = Cint(0)

# C-callable dummy callback
const dummy_callback_c = @cfunction(dummy_callback, Cint, (Ptr{p4est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t}))

#############################################################################
# Main program
#############################################################################

# Create a connectivity structure for the unit square.
unitsquare_connectivity = p4est_connectivity_new_unitsquare()
@test unitsquare_connectivity != C_NULL
@test p4est_connectivity_is_valid(unitsquare_connectivity) == 1
@test p4est_connectivity_complete(unitsquare_connectivity) == nothing
@test p4est_connectivity_reduce(unitsquare_connectivity) == nothing

# Test new connectivities
for test_conn in [p4est_connectivity_new_periodic(),
                p4est_connectivity_new_rotwrap(),
                p4est_connectivity_new_corner(),
                p4est_connectivity_new_pillow(),
                p4est_connectivity_new_moebius(),
                p4est_connectivity_new_star(),
                p4est_connectivity_new_cubed(),
                p4est_connectivity_new_brick(3, 18, 0, 1)]
    @test test_conn != C_NULL
    @test p4est_connectivity_is_valid(test_conn) == 1
    @test p4est_connectivity_is_equivalent(unitsquare_connectivity, test_conn) == 0
    @test p4est_connectivity_destroy(test_conn) == nothing
end

# Create a new forest
unitsquare_forest = p4est_new(mpicomm, unitsquare_connectivity, 0, C_NULL, C_NULL)
@test unitsquare_forest != C_NULL

# Create/Destroy geometry
unitsquare_geometry = p4est_geometry_new_connectivity(unitsquare_forest)
@test unitsquare_geometry != C_NULL
@test p4est_geometry_destroy(unitsquare_geometry) == nothing

# Copy/Destroy forest
unitsquare_forest_copy =  p4est_copy(unitsquare_forest, 1)
@test p4est_destroy(unitsquare_forest_copy) == nothing
unitsquare_forest_copy =  p4est_copy_ext(unitsquare_forest, 1, 0)
@test p4est_destroy(unitsquare_forest_copy) == nothing

# Operations on forests
@test p4est_partition(unitsquare_forest, 0, C_NULL) == nothing
@test p4est_partition_ext(unitsquare_forest, 0, C_NULL) == 0
@test p4est_balance(unitsquare_forest, P4est_wrapper.P4EST_CONNECT_FULL, C_NULL) == nothing
@test p4est_balance_ext(unitsquare_forest, P4est_wrapper.P4EST_CONNECT_FULL, C_NULL, C_NULL) == nothing
@test p4est_refine(unitsquare_forest, 1, dummy_callback_c, C_NULL) == nothing
@test p4est_refine_ext(unitsquare_forest, 1, -1, dummy_callback_c, C_NULL, C_NULL) == nothing
@test p4est_coarsen(unitsquare_forest, 1, dummy_callback_c, C_NULL) == nothing
@test p4est_coarsen_ext(unitsquare_forest, 1, 0, dummy_callback_c, C_NULL, C_NULL) == nothing
@test p4est_revision(unitsquare_forest) == 0
@test typeof(p4est_checksum(unitsquare_forest)) == UInt32
@test p4est_connect_type_int(P4est_wrapper.P4EST_CONNECT_FACE) == 1
@test unsafe_string(p4est_connect_type_string(P4est_wrapper.P4EST_CONNECT_FACE)) == "FACE"

# Destroy the forest
@test p4est_destroy(unitsquare_forest) == nothing
# Destroy the connectivity
@test p4est_connectivity_destroy(unitsquare_connectivity) == nothing

# Finalize MPI if initialized and session is not interactive
if (MPI.Initialized() && !isinteractive())
    MPI.Finalize()
end


text:
using LinearAlgebra: norm

"""
    euclidean_distance(a, b, [datum = wgs84])

The straight-line distance between points `a` and `b`. Uses `datum` to perform
transformations as necessary, and unlike other *Geodesy* functions, this
defaults to use the common WGS-84 datum for convenience.
"""
euclidean_distance(a::T, b::T) where {T <: AbstractVector} = norm(a-b)

# Automatic transformations to ECEF
# Uses wgs84 datum as a default. In most cases, the datum choice will only make
# a small difference to the answer. Nevertheless, is this acceptable?
euclidean_distance(a::T, b::T, datum) where {T <: AbstractVector} = euclidean_distance(a, b)
euclidean_distance(a::LLA, b, datum = wgs84) = euclidean_distance(ECEFfromLLA(datum)(a), b, datum)
euclidean_distance(a::ECEF, b::LLA, datum = wgs84) = euclidean_distance(a, ECEFfromLLA(datum)(b), datum)
euclidean_distance(a::LatLon, b, datum = wgs84) = euclidean_distance(ECEFfromLLA(datum)(LLA(a)), b, datum)
euclidean_distance(a::ECEF, b::LatLon, datum = wgs84) = euclidean_distance(a, ECEFfromLLA(datum)(LLA(b)), datum)
euclidean_distance(a::UTMZ, b, datum = wgs84) = euclidean_distance(ECEFfromUTMZ(datum)(a), b, datum)
euclidean_distance(a::ECEF, b::UTMZ, datum = wgs84) = euclidean_distance(a, ECEFfromUTMZ(datum)(b), datum)

"""
    euclidean_distance(utm1, utm2, zone, isnorth, [datum = wgs84])
    euclidean_distance(a, utm2, zone, isnorth, [datum = wgs84])
    euclidean_distance(utm1, b, zone, isnorth, [datum = wgs84])

If one or both points are UTM, we need the zone (and particularly the hemisphere,
isnorth = true/false) to determine the Euclidean distance.
"""
euclidean_distance(a::UTM, b::UTM, zone::Integer, isnorth::Bool, datum = wgs84) = euclidean_distance(ECEFfromUTM(zone, isnorth, datum)(a), ECEFfromUTM(zone, isnorth, datum)(b), datum)
euclidean_distance(a, b::UTM, zone::Integer, isnorth::Bool, datum = wgs84) = euclidean_distance(a, ECEFfromUTM(zone, isnorth, datum)(b), datum)
euclidean_distance(a::UTM, b, zone::Integer, isnorth::Bool, datum = wgs84) = euclidean_distance(ECEFfromUTM(zone, isnorth, datum)(a), b, datum)


text:
"""
    consecutive(k::Integer, sepp::SEPP)::AbstractVector{<:Real}

Compute `p_1,...,p_k` where `p_i` is the probability of an event on `i` knowing there has been events on `1,...,i-1`.
"""
function consecutive end

function consecutive(k::Integer, sepp::SEPPExpKern)::AbstractVector{<:Real}

    μ = sepp.μ
    ϕ = sepp.ϕ
    γ = sepp.γ

    times = 1:k
    ts = TimeSeries(times)

    vol = volfunc(times, ts, γ)
    intens = μ .+ ϕ .* vol
    prob = 1 .- exp.(-intens )

    return prob
end

#= TODO consecutive qui prend en compte une fonction impact pour les marques
function consecutive(k::Int, sempp::SEMPPExpKern)::AbstractVector{<:Real}

    μ = sempp.μ
    ϕ = sempp.ϕ
    γ = sempp.γ
    δ = sempp.δ
    markdens = sempp.markdens
    ξ = sempp.ξ
    α = sempp.α
    β = sempp.β
    κ = sempp.κ

    times = 1:k
    ts = TimeSeries(times)

    vol = volfunc(times, ts, γ)
    intens = μ .+ ϕ .* vol
    prob = 1 .- exp.(-intens )

    return prob
end
=#

"""
    markov_return_period(sepp::SEPP, r::Integer)::Real

Compute the expected return period of a run of r consecutive events for the model based on a markov chain.

See Li2021 4.2.
"""
function markov_return_period(sepp::SEPP, r::Integer)::Real # seulement pour SEPPExpKern pour l'instant
    p = consecutive(r+1,sepp)
    return (1 + sum(cumprod(p)[1:(r-1)]) + prod(p[1:r]) / (1 - p[r+1])) / (1 - sum((1 .- p[1:r]) .* cumprod(vcat([1], p[1:r-1]))))
end


text:
using Test, QuadGK, GrapheneQFT, LinearAlgebra

## Variables used for testing
a1 = GrapheneCoord(0, 17, A)
a2 = GrapheneCoord(1, 1, A)
a3 = GrapheneCoord(3, -2, B)
a4 = GrapheneCoord(2, 7, A)
a5 = GrapheneCoord(0, 1, A)

s1 = GrapheneState(a1, SpinUp)
s2 = GrapheneState(a2, SpinDown)
s3 = GrapheneState(a3, SpinUp)
s4 = GrapheneState(a4, SpinDown)
s5 = GrapheneState(a5, SpinUp)

ϵ1 = 1.7
ϵ2 = -0.2

V1 = 0.4
V2 = 2.0
V3 = -7.1

jx1 = 1.0
jy1 = 0.3
jz1 = 1.0

imp1 = ImpurityState(ϵ1, [(V1, a4), (V2, a3)])
imp2 = ImpurityState(ϵ2, [(V3, a5)])

sp1 = LocalSpin(jx1, jy1, jz1, a2)

c1 = 3.0 + 0im
c2 = 5.0 + 1im
c3 = -9.0 + 0im

h1 = Hopping(a1, a1, c1)
h2 = Hopping(a2, a1, c2)
h3 = Hopping(a3, a1, c3)
## Tests
@testset "Basic definitions" begin
    @test map(
        x -> GrapheneCoord(x.u, x.v, !x.sublattice),
        [a1, a2, a3, a4, a5],
    ) == [
        GrapheneCoord(0, 17, B),
        GrapheneCoord(1, 1, B),
        GrapheneCoord(3, -2, A),
        GrapheneCoord(2, 7, B),
        GrapheneCoord(0, 1, B),
    ]

    @test map(x -> GrapheneState(x.coord, !x.spin), [s1, s2, s3, s4, s5]) == [
        GrapheneState(a1, SpinDown),
        GrapheneState(a2, SpinUp),
        GrapheneState(a3, SpinDown),
        GrapheneState(a4, SpinUp),
        GrapheneState(a5, SpinDown),
    ]

    @test isless(a2, a1) == true
    @test isless(a3, a1) == false
    @test isless(a2, a5) == false
    @test isless(a1, a5) == false

    @test isless(s2, s1) == true
    @test isless(s3, s1) == false
    @test isless(s2, s5) == true
    @test isless(s1, s5) == false

    @test repr(a1) == "|0, 17, A⟩"
    @test repr(s3) == "|3, -2, B⟩⊗|SpinUp⟩"

end

@testset "Helper functions" begin
    @test sort(graphene_neighbors(a4)) == [
        GrapheneCoord(2, 7, B),
        GrapheneCoord(3, 7, B),
        GrapheneCoord(2, 8, B),
    ]

    @test sort(graphene_neighbors(a3)) == [
        GrapheneCoord(3, -3, A),
        GrapheneCoord(2, -2, A),
        GrapheneCoord(3, -2, A),
    ]
    #
    # @test vcat(a2, sort(graphene_neighbors(a2))) ==
    #       sort(graphene_multiple_neighbors(a2, 1))

    @test (
        crystal_to_cartesian(a1) .≈ [
            GrapheneQFT.graphene_d1[1] * 0 + GrapheneQFT.graphene_d2[1] * 17,
            GrapheneQFT.graphene_d1[2] * 0 + GrapheneQFT.graphene_d2[2] * 17,
        ]
    ) |> prod
    @test (
        crystal_to_cartesian(a3) .≈ [
            GrapheneQFT.graphene_d1[1] * 3 + GrapheneQFT.graphene_d2[1] * (-2),
            GrapheneQFT.graphene_d1[2] * 3 +
            GrapheneQFT.graphene_d2[2] * (-2) +
            GrapheneQFT.sublattice_shift,
        ]
    ) |> prod
end

@testset "No Impurity States" begin
    rand_num = ((rand() - 1 / 2) + 1im * (rand() - 1 / 2)) * 5
    my_system = mkGrapheneSystem(0.0, 0.0, [sp1, h1, h2, h3])
    @test my_system.T == 0.0
    @test my_system.μ == 0.0
    @test my_system.imps == Vector{Float64}[]
    @test my_system.scattering_states == vcat(
        [GrapheneState(x, SpinUp) for x in sort([a1, a2, a3])],
        [GrapheneState(x, SpinDown) for x in sort([a1, a2, a3])],
    )
    @test my_system.V == Array{Float64,2}(undef, 0, 0)
    @test my_system.Δ == [
        jz1 c2 0 (jx1-1im*jy1) 0 0
        conj(c2) c1 conj(c3) 0 0 0
        0 c3 0 0 0 0
        (jx1+1im*jy1) 0 0 (-jz1) c2 0
        0 0 0 conj(c2) c1 conj(c3)
        0 0 0 0 c3 0
    ]

    scattering_pair = [(s4, s5)]
    scattering_pair_rev = [(s5, s4)]

    @test G_R(rand_num, scattering_pair, my_system) |> real ≈
          G_R(conj.(rand_num), scattering_pair_rev, my_system) |> real

    @test G_R(rand_num, scattering_pair, my_system) |> imag ≈
          -G_R(conj.(rand_num), scattering_pair_rev, my_system) |> imag

    @test δG_R(rand_num, [(s1, s2), (s3, s4)], my_system) == G_R(rand_num, [(s1, s2), (s3, s4)], my_system)

end

@testset "Only Impurities" begin
    rand_num = ((rand() - 1 / 2) + 1im * (rand() - 1 / 2)) * 5
    my_system = mkGrapheneSystem(0.0, 0.0, Vector{Defect}([imp1, imp2]))
    @test my_system.T == 0.0
    @test my_system.μ == 0.0
    @test my_system.imps == [ϵ1, ϵ2]
    @test my_system.scattering_states == vcat(
        [GrapheneState(x, SpinUp) for x in sort([a3, a4, a5])],
        [GrapheneState(x, SpinDown) for x in sort([a3, a4, a5])],
    )
    @test my_system.V == [
        0 V3 0 0
        V1 0 0 0
        V2 0 0 0
        0 0 0 V3
        0 0 V1 0
        0 0 V2 0
    ]

    @test my_system.Δ == [
        0 0 0 0 0 0
        0 0 0 0 0 0
        0 0 0 0 0 0
        0 0 0 0 0 0
        0 0 0 0 0 0
        0 0 0 0 0 0
    ]

    scattering_pair = [(s1, s5)]
    scattering_pair_rev = [(s5, s1)]

    @test G_R(rand_num, scattering_pair, my_system) |> real ≈
          G_R(conj.(rand_num), scattering_pair_rev, my_system) |> real

    @test G_R(rand_num, scattering_pair, my_system) |> imag ≈
          -G_R(conj.(rand_num), scattering_pair_rev, my_system) |> imag

    @test G_R(rand_num, [(s1, s2)], my_system)[1] ≈ 0

    @test δΓ(rand_num, my_system) |> Diagonal ≈ conj.(δΓ(conj.(rand_num), my_system) |> Diagonal)

    @test Γ(rand_num, my_system) |> Diagonal ≈ conj.(Γ(conj.(rand_num), my_system) |> Diagonal)

    @test δρ_R_graphene(GrapheneState(a1, SpinUp), my_system) == δρ_R_graphene(GrapheneState(a1, SpinDown), my_system)

end

@testset "Impurities and Perturbation" begin
    rand_num = ((rand() - 1 / 2) + 1im * (rand() - 1 / 2)) * 5
    my_system = mkGrapheneSystem(0.0, 0.0, [sp1, h1, h2, h3, imp1, imp2])
    @test my_system.T == 0.0
    @test my_system.μ == 0.0
    @test my_system.imps == [ϵ1, ϵ2]
    @test my_system.scattering_states == vcat(
        [GrapheneState(x, SpinUp) for x in sort([a1, a2, a3, a4, a5])],
        [GrapheneState(x, SpinDown) for x in sort([a1, a2, a3, a4, a5])],
    )

    @test my_system.V == [
        0 V3 0 0
        0 0 0 0
        V1 0 0 0
        0 0 0 0
        V2 0 0 0
        0 0 0 V3
        0 0 0 0
        0 0 V1 0
        0 0 0 0
        0 0 V2 0
    ]

    @test my_system.Δ == [
        0 0 0 0 0 0 0 0 0 0
        0 jz1 0 c2 0 0 (jx1-1im*jy1) 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 conj(c2) 0 c1 conj(c3) 0 0 0 0 0
        0 0 0 c3 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
        0 (jx1+1im*jy1) 0 0 0 0 (-jz1) 0 c2 0
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 conj(c2) 0 c1 conj(c3)
        0 0 0 0 0 0 0 0 c3 0
    ]

    scattering_pair = [(s1, s5), (s4, s2)]
    scattering_pair_rev = [(s5, s1), (s2, s4)]

    @test G_R(rand_num, scattering_pair, my_system) |> real ≈
          G_R(conj.(rand_num), scattering_pair_rev, my_system) |> real

    @test G_R(rand_num, scattering_pair, my_system) |> imag ≈
          -G_R(conj.(rand_num), scattering_pair_rev, my_system) |> imag

    @test δΓ(rand_num, my_system) |> Diagonal ≈ conj.(δΓ(conj.(rand_num), my_system) |> Diagonal)

end

@testset "No defects present" begin
    rand_num = ((rand() - 1 / 2) + 1im * (rand() - 1 / 2)) * 5
    my_system = mkGrapheneSystem(0.0, 0.0, Vector{Defect}([]))

    @test my_system.Δ == Array{GrapheneCoord}(undef, 0, 0)
    scattering_pair = [(s1, s5), (s4, s2)]
    @test δG_R(rand_num, scattering_pair, my_system) ==
          [(0.0 + 0.0im), (0.0 + 0.0im)]
    @test δρ_R_graphene(s2, my_system) == 0.0
    @test δρ_R_graphene(s3, my_system) == 0.0
end


text:
export isiban

function isiban(str::AbstractString)::Bool
    ibanReg = Dict([
        ("AD", r"^(AD[0-9]{2})\d{8}[A-Z0-9]{12}$")
        ("AE", r"^(AE[0-9]{2})\d{3}\d{16}$")
        ("AL", r"^(AL[0-9]{2})\d{8}[A-Z0-9]{16}$")
        ("AT", r"^(AT[0-9]{2})\d{16}$")
        ("AZ", r"^(AZ[0-9]{2})[A-Z0-9]{4}\d{20}$")
        ("BA", r"^(BA[0-9]{2})\d{16}$")
        ("BE", r"^(BE[0-9]{2})\d{12}$")
        ("BG", r"^(BG[0-9]{2})[A-Z]{4}\d{6}[A-Z0-9]{8}$")
        ("BH", r"^(BH[0-9]{2})[A-Z]{4}[A-Z0-9]{14}$")
        ("BR", r"^(BR[0-9]{2})\d{23}[A-Z]{1}[A-Z0-9]{1}$")
        ("BY", r"^(BY[0-9]{2})[A-Z0-9]{4}\d{20}$")
        ("CH", r"^(CH[0-9]{2})\d{5}[A-Z0-9]{12}$")
        ("CR", r"^(CR[0-9]{2})\d{18}$")
        ("CY", r"^(CY[0-9]{2})\d{8}[A-Z0-9]{16}$")
        ("CZ", r"^(CZ[0-9]{2})\d{20}$")
        ("DE", r"^(DE[0-9]{2})\d{18}$")
        ("DK", r"^(DK[0-9]{2})\d{14}$")
        ("DO", r"^(DO[0-9]{2})[A-Z]{4}\d{20}$")
        ("EE", r"^(EE[0-9]{2})\d{16}$")
        ("EG", r"^(EG[0-9]{2})\d{25}$")
        ("ES", r"^(ES[0-9]{2})\d{20}$")
        ("FI", r"^(FI[0-9]{2})\d{14}$")
        ("FO", r"^(FO[0-9]{2})\d{14}$")
        ("FR", r"^(FR[0-9]{2})\d{10}[A-Z0-9]{11}\d{2}$")
        ("GB", r"^(GB[0-9]{2})[A-Z]{4}\d{14}$")
        ("GE", r"^(GE[0-9]{2})[A-Z0-9]{2}\d{16}$")
        ("GI", r"^(GI[0-9]{2})[A-Z]{4}[A-Z0-9]{15}$")
        ("GL", r"^(GL[0-9]{2})\d{14}$")
        ("GR", r"^(GR[0-9]{2})\d{7}[A-Z0-9]{16}$")
        ("GT", r"^(GT[0-9]{2})[A-Z0-9]{4}[A-Z0-9]{20}$")
        ("HR", r"^(HR[0-9]{2})\d{17}$")
        ("HU", r"^(HU[0-9]{2})\d{24}$")
        ("IE", r"^(IE[0-9]{2})[A-Z0-9]{4}\d{14}$")
        ("IL", r"^(IL[0-9]{2})\d{19}$")
        ("IQ", r"^(IQ[0-9]{2})[A-Z]{4}\d{15}$")
        ("IR", r"^(IR[0-9]{2})0\d{2}0\d{18}$")
        ("IS", r"^(IS[0-9]{2})\d{22}$")
        ("IT", r"^(IT[0-9]{2})[A-Z]{1}\d{10}[A-Z0-9]{12}$")
        ("JO", r"^(JO[0-9]{2})[A-Z]{4}\d{22}$")
        ("KW", r"^(KW[0-9]{2})[A-Z]{4}[A-Z0-9]{22}$")
        ("KZ", r"^(KZ[0-9]{2})\d{3}[A-Z0-9]{13}$")
        ("LB", r"^(LB[0-9]{2})\d{4}[A-Z0-9]{20}$")
        ("LC", r"^(LC[0-9]{2})[A-Z]{4}[A-Z0-9]{24}$")
        ("LI", r"^(LI[0-9]{2})\d{5}[A-Z0-9]{12}$")
        ("LT", r"^(LT[0-9]{2})\d{16}$")
        ("LU", r"^(LU[0-9]{2})\d{3}[A-Z0-9]{13}$")
        ("LV", r"^(LV[0-9]{2})[A-Z]{4}[A-Z0-9]{13}$")
        ("MC", r"^(MC[0-9]{2})\d{10}[A-Z0-9]{11}\d{2}$")
        ("MD", r"^(MD[0-9]{2})[A-Z0-9]{20}$")
        ("ME", r"^(ME[0-9]{2})\d{18}$")
        ("MK", r"^(MK[0-9]{2})\d{3}[A-Z0-9]{10}\d{2}$")
        ("MR", r"^(MR[0-9]{2})\d{23}$")
        ("MT", r"^(MT[0-9]{2})[A-Z]{4}\d{5}[A-Z0-9]{18}$")
        ("MU", r"^(MU[0-9]{2})[A-Z]{4}\d{19}[A-Z]{3}$")
        ("MZ", r"^(MZ[0-9]{2})\d{21}$")
        ("NL", r"^(NL[0-9]{2})[A-Z]{4}\d{10}$")
        ("NO", r"^(NO[0-9]{2})\d{11}$")
        ("PK", r"^(PK[0-9]{2})[A-Z0-9]{4}\d{16}$")
        ("PL", r"^(PL[0-9]{2})\d{24}$")
        ("PS", r"^(PS[0-9]{2})[A-Z0-9]{4}\d{21}$")
        ("PT", r"^(PT[0-9]{2})\d{21}$")
        ("QA", r"^(QA[0-9]{2})[A-Z]{4}[A-Z0-9]{21}$")
        ("RO", r"^(RO[0-9]{2})[A-Z]{4}[A-Z0-9]{16}$")
        ("RS", r"^(RS[0-9]{2})\d{18}$")
        ("SA", r"^(SA[0-9]{2})\d{2}[A-Z0-9]{18}$")
        ("SC", r"^(SC[0-9]{2})[A-Z]{4}\d{20}[A-Z]{3}$")
        ("SE", r"^(SE[0-9]{2})\d{20}$")
        ("SI", r"^(SI[0-9]{2})\d{15}$")
        ("SK", r"^(SK[0-9]{2})\d{20}$")
        ("SM", r"^(SM[0-9]{2})[A-Z]{1}\d{10}[A-Z0-9]{12}$")
        ("SV", r"^(SV[0-9]{2})[A-Z0-9]{4}\d{20}$")
        ("TL", r"^(TL[0-9]{2})\d{19}$")
        ("TN", r"^(TN[0-9]{2})\d{20}$")
        ("TR", r"^(TR[0-9]{2})\d{5}[A-Z0-9]{17}$")
        ("UA", r"^(UA[0-9]{2})\d{6}[A-Z0-9]{19}$")
        ("VA", r"^(VA[0-9]{2})\d{18}$")
        ("VG", r"^(VG[0-9]{2})[A-Z0-9]{4}\d{16}$")
        ("XK", r"^(XK[0-9]{2})\d{16}$")
    ])

    # Check format

    stripped = replace(str, r"[\s\-]+"i => "")
    countryCode = uppercase(stripped[1:2])
    if !(countryCode in keys(ibanReg))
        return false
    end

    if !contains(stripped, ibanReg[countryCode])
        return false
    end

    # Check checksum
    stripped = uppercase(replace(str, r"[^A-Z0-9]+"i => ""))
    rearranged = str[5:end] * str[1:4]
    alphaCapsReplacedWithDigits = join(map(((index, value), ) -> contains(string(value), r"[A-Z]") ? Int(codepoint(value)) - 55 : value, enumerate(rearranged)))
    splitted = SubString.(alphaCapsReplacedWithDigits, findall(r"\d{1,7}", alphaCapsReplacedWithDigits))
    remainder = parse(Int, reduce((acc, val) -> string(parse(Int, acc * val) % 97), splitted))

    return remainder === 1
end


text:
using Peacock, PyPlot, LinearAlgebra
a = 1; a1 = [1, 0]; a2 = [0, 1]
w = 0.17 # 柱子边长
d = w/2 # w/2 中间， a/2 - w/2 四周 a/4 中心

eps0 = 1; mu0 = 1
mu1 = [14 12.4im 0; -12.4im 14 0; 0 0 1]
eps1 = 15

function epf(x, y)
    if abs(x - d) < w/2 && abs(y - d) < w/2 || abs(x + d) < w/2 && abs(y - d) < w/2 || abs(x - d) < w/2 && abs(y + d) < w/2 || abs(x + d) < w/2 && abs(y + d) < w/2
        return eps1
    else
        return eps0
    end
end

function muf(x, y)
    if abs(x - d) < w/2 && abs(y - d) < w/2 || abs(x + d) < w/2 && abs(y - d) < w/2 || abs(x - d) < w/2 && abs(y + d) < w/2 || abs(x + d) < w/2 && abs(y + d) < w/2
        return mu1
    else
        return mu0
    end
end

#epf(x, y) = x^2 + y^2 < 0.11^2 ? 15 : 1
#muf(x, y) = x^2 + y^2 < 0.11^2 ? mu1 : 1

d1 = d2 = 0.01
geometry = Geometry(epf, muf, a1, a2, d1, d2, TM)
ks = [
    BrillouinZoneCoordinate(-0.5, -0.5, "-X"),
    BrillouinZoneCoordinate(-0.5, 0.0, "Γ"),
    BrillouinZoneCoordinate(-0.5, 0.5, "X")
]
solver = Peacock.FDFD.Solver(geometry, 2*[d1, d2])

dk_inner = nothing
dk_outer = 0.25
bands = 1:2
figure()
plot_wilson_loop_winding(solver, ks, TM, bands, dk_outer=dk_outer, dk_inner=dk_inner,delta_brillouin_zone=BrillouinZoneCoordinate(1,0))


text:
include("./deeper-cut/mpii.annotation.reader.jl")
using Plots; default(fmt = :png)

function draw_plots(results)
    trnloss,tstloss = Array{Float32}(results[2,:]), Array{Float32}(results[3,:])
    plot_loss = plot([trnloss,tstloss],ylim=(.0,.4),labels=["trnloss" "valloss"],xlabel="Epochs",ylabel="Loss")

        trnacc,tstacc = Array{Float32}(results[6,:]), Array{Float32}(results[7,:])
    plot_acc = plot([trnacc,tstacc],ylim=(.0,1.0),labels=["tracc" "valacc"],xlabel="Epochs",ylabel="Accuracy")

    trnerr,tsterr = Array{Float32}(results[4,:]), Array{Float32}(results[5,:])
    plot_error = plot([trnerr,tsterr],ylim=(.0,.12),labels=["trnerr" "valerr"],xlabel="Epochs",ylabel="Error")

    trnacc_pckh, tstacc_pckh = Array{Float32}(results[8,:]), Array{Float32}(results[9,:])
    plot_acc_pckh = plot([trnacc_pckh, tstacc_pckh],ylim=(.0,1.0),labels=["tracc_pckh" "valacc_pckh"],xlabel="Epochs",ylabel="Accuracy PCKh@$(PCKh_range)")

    foreach(display,
            [plot_loss,
            plot_acc,
            plot_acc_pckh,
            plot_error])
end


text:
include("utils.jl")

@testset "SkewSymmetricMatrices" begin
    M = SkewSymmetricMatrices(3, ℝ)
    A = [1 2 3; 4 5 6; 7 8 9]
    A_skewsym = [0 -2 -3; 2 0 1; 3 -1 0]
    A_skewsym2 = [0 -2 -3; 2 0 1; 3 -1 0]
    B_skewsym = [0 -2 -3; 2 0 -1; 3 1 0]
    M_complex = SkewSymmetricMatrices(3, ℂ)
    A_skewsym_complex = [1.0im -2.0 -3.0; 2.0 0.0 1.0; 3.0 -1.0 0.0]
    B_skewsym_complex = [2.0im -2.0 - 2.0im -3.0; 2.0 - 2.0im 0.0 -1.0; 3.0 1.0 0.0]
    @test repr(M_complex) == "SkewSymmetricMatrices(3, ℂ)"
    C = [0 -1 im; 1 0 -im; im -im 0]
    D = [1 0; 0 1]
    X = zeros(3, 3)
    @testset "Real Skew-Symmetric Matrices Basics" begin
        @test repr(M) == "SkewSymmetricMatrices(3, ℝ)"
        @test representation_size(M) == (3, 3)
        @test base_manifold(M) === M
        @test typeof(get_embedding(M)) === Euclidean{Tuple{3,3},ℝ}
        @test check_manifold_point(M, B_skewsym) === nothing
        @test_throws DomainError is_manifold_point(M, A, true)
        @test_throws DomainError is_manifold_point(M, C, true)
        @test_throws DomainError is_manifold_point(M, D, true)
        @test check_tangent_vector(M, B_skewsym, B_skewsym) === nothing
        @test_throws DomainError is_tangent_vector(M, B_skewsym, A, true)
        @test_throws DomainError is_tangent_vector(M, A, B_skewsym, true)
        @test_throws DomainError is_tangent_vector(M, B_skewsym, D, true)
        @test_throws DomainError is_tangent_vector(
            M,
            B_skewsym,
            1 * im * zero_tangent_vector(M, B_skewsym),
            true,
        )
        @test manifold_dimension(M) == 3
        @test manifold_dimension(M_complex) == 9
        @test A_skewsym2 == project!(M, A_skewsym, A_skewsym)
        @test A_skewsym2 == project(M, A_skewsym, A_skewsym)
        A_sym3 = similar(A_skewsym)
        embed!(M, A_sym3, A_skewsym)
        A_sym4 = embed(M, A_skewsym)
        @test A_sym3 == A_skewsym
        @test A_sym4 == A_skewsym
    end
    types = [Matrix{Float64}, ]
    TEST_FLOAT32 && push!(types, Matrix{Float32})
    TEST_STATIC_SIZED && push!(types, MMatrix{3, 3, Float64})
    bases = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))
    for T in types
        pts = [convert(T, A_skewsym), convert(T, B_skewsym), convert(T, X)]
        @testset "Type $T" begin
            test_manifold(
                M,
                pts,
                test_injectivity_radius = false,
                test_reverse_diff = isa(T, Vector),
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true,
                basis_types_vecs = (
                    DiagonalizingOrthonormalBasis(log(M, pts[1], pts[2])),
                    bases...,
                ),
                basis_types_to_from = bases,
            )
        end
    end
    complex_types = [ Matrix{ComplexF64}, ]
    TEST_FLOAT32 && push!(complex_types, Matrix{ComplexF32})
    TEST_STATIC_SIZED && push!(compley_types, MMatrix{3, 3, ComplexF64})
    for T in complex_types
        pts_complex =
            [convert(T, A_skewsym_complex), convert(T, B_skewsym_complex), convert(T, X)]
        @testset "Type $T" begin
            test_manifold(
                M_complex,
                pts_complex,
                test_injectivity_radius = false,
                test_reverse_diff = isa(T, Vector),
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true,
                basis_types_vecs = (DefaultOrthonormalBasis(),),
                basis_types_to_from = (DefaultOrthonormalBasis(),),
            )
            @test isapprox(
                -pts_complex[1],
                exp(M, pts_complex[1], log(M, pts_complex[1], -pts_complex[1])),
            )
        end # testset type $T
    end # for
end # test SymmetricMatrices