text:
using Mocking

check_dnsavailability_1_patch = @patch function check_dnsavailability(a...; b...)
    return Dict{String, Any}("Available" => true, "FullyQualifiedCNAME" => "my-cname.us-west-2.elasticbeanstalk.com")
end

create_application_1_patch = @patch function create_application(a...; b...)
    return Dict{String, Any}("Application" => Dict{String, Any}("DateCreated" => "2015-02-12T18:32:21.181Z", "DateUpdated" => "2015-02-12T18:32:21.181Z", "Description" => "my application", "ApplicationName" => "my-app", "ConfigurationTemplates" => Any[]))
end

create_application_version_1_patch = @patch function create_application_version(a...; b...)
    return Dict{String, Any}("ApplicationVersion" => Dict{String, Any}("SourceBundle" => Dict{String, Any}("S3Key" => "sample.war", "S3Bucket" => "my-bucket"), "DateCreated" => "2015-02-03T23:01:25.412Z", "DateUpdated" => "2015-02-03T23:01:25.412Z", "Description" => "my-app-v1", "VersionLabel" => "v1", "ApplicationName" => "my-app"))
end

create_configuration_template_1_patch = @patch function create_configuration_template(a...; b...)
    return Dict{String, Any}("TemplateName" => "my-app-v1", "DateCreated" => "2015-08-12T18:40:39Z", "DateUpdated" => "2015-08-12T18:40:39Z", "ApplicationName" => "my-app", "SolutionStackName" => "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 8 Java 8")
end

create_environment_1_patch = @patch function create_environment(a...; b...)
    return Dict{String, Any}("DateCreated" => "2015-02-03T23:04:54.479Z", "VersionLabel" => "v1", "SolutionStackName" => "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 8 Java 8", "DateUpdated" => "2015-02-03T23:04:54.479Z", "Status" => "Launching", "Tier" => Dict{String, Any}("Type" => "Standard", "Version" => " ", "Name" => "WebServer"), "EnvironmentName" => "my-env", "CNAME" => "my-app.elasticbeanstalk.com", "Health" => "Grey", "ApplicationName" => "my-app", "EnvironmentId" => "e-izqpassy4h")
end

create_storage_location_1_patch = @patch function create_storage_location(a...; b...)
    return Dict{String, Any}("S3Bucket" => "elasticbeanstalk-us-west-2-0123456789012")
end

describe_application_versions_1_patch = @patch function describe_application_versions(a...; b...)
    return Dict{String, Any}("ApplicationVersions" => Any[Dict{String, Any}("SourceBundle" => Dict{String, Any}("S3Key" => "my-app/5026-stage-150723_224258.war", "S3Bucket" => "elasticbeanstalk-us-west-2-015321684451"), "DateCreated" => "2015-07-23T01:32:26.079Z", "DateUpdated" => "2015-07-23T01:32:26.079Z", "Description" => "update cover page", "VersionLabel" => "v2", "ApplicationName" => "my-app"), Dict{String, Any}("SourceBundle" => Dict{String, Any}("S3Key" => "my-app/5026-stage-150723_222618.war", "S3Bucket" => "elasticbeanstalk-us-west-2-015321684451"), "DateCreated" => "2015-07-23T22:26:10.816Z", "DateUpdated" => "2015-07-23T22:26:10.816Z", "Description" => "initial version", "VersionLabel" => "v1", "ApplicationName" => "my-app")])
end

describe_applications_1_patch = @patch function describe_applications(a...; b...)
    return Dict{String, Any}("Applications" => Any[Dict{String, Any}("Versions" => Any["Sample Application"], "DateCreated" => "2015-08-13T21:05:44.376Z", "DateUpdated" => "2015-08-13T21:05:44.376Z", "ApplicationName" => "ruby", "ConfigurationTemplates" => Any[]), Dict{String, Any}("Versions" => Any["Sample Application"], "DateCreated" => "2015-08-13T19:05:43.637Z", "DateUpdated" => "2015-08-13T19:05:43.637Z", "Description" => "Application created from the EB CLI using \"eb init\"", "ApplicationName" => "pythonsample", "ConfigurationTemplates" => Any[]), Dict{String, Any}("Versions" => Any["add elasticache", "First Release"], "DateCreated" => "2015-08-06T17:50:02.486Z", "DateUpdated" => "2015-08-06T17:50:02.486Z", "ApplicationName" => "nodejs-example", "ConfigurationTemplates" => Any[])])
end

describe_configuration_options_1_patch = @patch function describe_configuration_options(a...; b...)
    return Dict{String, Any}("Options" => Any[Dict{String, Any}("DefaultValue" => "30", "MaxValue" => 300, "ValueType" => "Scalar", "MinValue" => 5, "Namespace" => "aws:elb:healthcheck", "Name" => "Interval", "ChangeSeverity" => "NoInterruption", "UserDefined" => false), Dict{String, Any}("DefaultValue" => "2000000", "ValueType" => "Scalar", "MinValue" => 0, "Namespace" => "aws:autoscaling:trigger", "Name" => "LowerThreshold", "ChangeSeverity" => "NoInterruption", "UserDefined" => false)])
end

describe_configuration_settings_1_patch = @patch function describe_configuration_settings(a...; b...)
    return Dict{String, Any}("ConfigurationSettings" => Any[Dict{String, Any}("OptionSettings" => Any[Dict{String, Any}("Value" => "Any", "ResourceName" => "AWSEBAutoScalingGroup", "OptionName" => "Availability Zones", "Namespace" => "aws:autoscaling:asg"), Dict{String, Any}("Value" => "360", "ResourceName" => "AWSEBAutoScalingGroup", "OptionName" => "Cooldown", "Namespace" => "aws:autoscaling:asg"), Dict{String, Any}("Value" => "20", "ResourceName" => "AWSEBLoadBalancer", "OptionName" => "ConnectionDrainingTimeout", "Namespace" => "aws:elb:policies"), Dict{String, Any}("Value" => "60", "ResourceName" => "AWSEBLoadBalancer", "OptionName" => "ConnectionSettingIdleTimeout", "Namespace" => "aws:elb:policies")], "DateCreated" => "2015-08-13T19:16:25Z", "DateUpdated" => "2015-08-13T23:30:07Z", "Description" => "Environment created from the EB CLI using \"eb create\"", "DeploymentStatus" => "deployed", "SolutionStackName" => "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 8 Java 8", "ApplicationName" => "my-app", "EnvironmentName" => "my-env")])
end

describe_environment_health_1_patch = @patch function describe_environment_health(a...; b...)
    return Dict{String, Any}("EnvironmentName" => "my-env", "InstancesHealth" => Dict{String, Any}("Pending" => 0, "NoData" => 0, "Degraded" => 0, "Severe" => 0, "Ok" => 1, "Info" => 0, "Unknown" => 0, "Warning" => 0), "HealthStatus" => "Ok", "ApplicationMetrics" => Dict{String, Any}("Latency" => Dict{String, Any}("P90" => 0.003, "P50" => 0.001, "P95" => 0.004, "P75" => 0.002, "P85" => 0.003, "P99" => 0.004, "P999" => 0.004, "P10" => 0.001), "RequestCount" => 45, "Duration" => 10, "StatusCodes" => Dict{String, Any}("Status3xx" => 0, "Status4xx" => 0, "Status2xx" => 45, "Status5xx" => 0)), "RefreshedAt" => "2015-08-20T21:09:18Z", "Causes" => Any[], "Color" => "Green")
end

describe_environment_resources_1_patch = @patch function describe_environment_resources(a...; b...)
    return Dict{String, Any}("EnvironmentResources" => Dict{String, Any}("AutoScalingGroups" => Any[Dict{String, Any}("Name" => "awseb-e-qu3fyyjyjs-stack-AWSEBAutoScalingGroup-QSB2ZO88SXZT")], "Instances" => Any[Dict{String, Any}("Id" => "i-0c91c786")], "LaunchConfigurations" => Any[Dict{String, Any}("Name" => "awseb-e-qu3fyyjyjs-stack-AWSEBAutoScalingLaunchConfiguration-1UUVQIBC96TQ2")], "Queues" => Any[], "Triggers" => Any[], "LoadBalancers" => Any[Dict{String, Any}("Name" => "awseb-e-q-AWSEBLoa-1EEPZ0K98BIF0")], "EnvironmentName" => "my-env"))
end

describe_environments_1_patch = @patch function describe_environments(a...; b...)
    return Dict{String, Any}("Environments" => Any[Dict{String, Any}("DateCreated" => "2015-08-07T20:48:49.599Z", "AbortableOperationInProgress" => false, "VersionLabel" => "7f58-stage-150812_025409", "SolutionStackName" => "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 8 Java 8", "EndpointURL" => "awseb-e-w-AWSEBLoa-1483140XB0Q4L-109QXY8121.us-west-2.elb.amazonaws.com", "DateUpdated" => "2015-08-12T18:16:55.019Z", "Status" => "Ready", "Tier" => Dict{String, Any}("Type" => "Standard", "Version" => " ", "Name" => "WebServer"), "EnvironmentName" => "my-env", "CNAME" => "my-env.elasticbeanstalk.com", "Health" => "Green", "ApplicationName" => "my-app", "EnvironmentId" => "e-rpqsewtp2j")])
end

describe_events_1_patch = @patch function describe_events(a...; b...)
    return Dict{String, Any}("Events" => Any[Dict{String, Any}("EventDate" => "2015-08-20T07:06:53.535Z", "Message" => "Environment health has transitioned from Info to Ok.", "ApplicationName" => "my-app", "Severity" => "INFO", "EnvironmentName" => "my-env"), Dict{String, Any}("EventDate" => "2015-08-20T07:06:02.049Z", "Message" => "Environment update completed successfully.", "ApplicationName" => "my-app", "Severity" => "INFO", "RequestId" => "b7f3960b-4709-11e5-ba1e-07e16200da41", "EnvironmentName" => "my-env"), Dict{String, Any}("EventDate" => "2015-08-13T19:16:27.561Z", "Message" => "Using elasticbeanstalk-us-west-2-012445113685 as Amazon S3 storage bucket for environment data.", "ApplicationName" => "my-app", "Severity" => "INFO", "RequestId" => "ca8dfbf6-41ef-11e5-988b-651aa638f46b", "EnvironmentName" => "my-env"), Dict{String, Any}("EventDate" => "2015-08-13T19:16:26.581Z", "Message" => "createEnvironment is starting.", "ApplicationName" => "my-app", "Severity" => "INFO", "RequestId" => "cdfba8f6-41ef-11e5-988b-65638f41aa6b", "EnvironmentName" => "my-env")])
end

describe_instances_health_1_patch = @patch function describe_instances_health(a...; b...)
    return Dict{String, Any}("InstanceHealthList" => Any[Dict{String, Any}("System" => Dict{String, Any}("LoadAverage" => Any[0, 0.02, 0.05], "CPUUtilization" => Dict{String, Any}("IOWait" => 0.2, "System" => 0.3, "User" => 1.5, "Nice" => 0.1, "Idle" => 97.8, "IRQ" => 0, "SoftIRQ" => 0.1)), "HealthStatus" => "Ok", "InstanceId" => "i-08691cc7", "LaunchedAt" => "2015-08-13T19:17:09Z", "ApplicationMetrics" => Dict{String, Any}("Latency" => Dict{String, Any}("P90" => 0.004, "P50" => 0.001, "P95" => 0.005, "P75" => 0.002, "P85" => 0.003, "P99" => 0.006, "P999" => 0.006, "P10" => 0), "RequestCount" => 48, "Duration" => 10, "StatusCodes" => Dict{String, Any}("Status3xx" => 0, "Status4xx" => 1, "Status2xx" => 47, "Status5xx" => 0)), "Causes" => Any[], "Color" => "Green")], "RefreshedAt" => "2015-08-20T21:09:08Z")
end

list_available_solution_stacks_1_patch = @patch function list_available_solution_stacks(a...; b...)
    return Dict{String, Any}("SolutionStackDetails" => Any[Dict{String, Any}("PermittedFileTypes" => Any["zip"], "SolutionStackName" => "64bit Amazon Linux 2015.03 v2.0.0 running Node.js")], "SolutionStacks" => Any["64bit Amazon Linux 2015.03 v2.0.0 running Node.js", "64bit Amazon Linux 2015.03 v2.0.0 running PHP 5.6", "64bit Amazon Linux 2015.03 v2.0.0 running PHP 5.5", "64bit Amazon Linux 2015.03 v2.0.0 running PHP 5.4", "64bit Amazon Linux 2015.03 v2.0.0 running Python 3.4", "64bit Amazon Linux 2015.03 v2.0.0 running Python 2.7", "64bit Amazon Linux 2015.03 v2.0.0 running Python", "64bit Amazon Linux 2015.03 v2.0.0 running Ruby 2.2 (Puma)", "64bit Amazon Linux 2015.03 v2.0.0 running Ruby 2.2 (Passenger Standalone)", "64bit Amazon Linux 2015.03 v2.0.0 running Ruby 2.1 (Puma)", "64bit Amazon Linux 2015.03 v2.0.0 running Ruby 2.1 (Passenger Standalone)", "64bit Amazon Linux 2015.03 v2.0.0 running Ruby 2.0 (Puma)", "64bit Amazon Linux 2015.03 v2.0.0 running Ruby 2.0 (Passenger Standalone)", "64bit Amazon Linux 2015.03 v2.0.0 running Ruby 1.9.3", "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 8 Java 8", "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 7 Java 7", "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 7 Java 6", "64bit Windows Server Core 2012 R2 running IIS 8.5", "64bit Windows Server 2012 R2 running IIS 8.5", "64bit Windows Server 2012 running IIS 8", "64bit Windows Server 2008 R2 running IIS 7.5", "64bit Amazon Linux 2015.03 v2.0.0 running Docker 1.6.2", "64bit Amazon Linux 2015.03 v2.0.0 running Multi-container Docker 1.6.2 (Generic)", "64bit Debian jessie v2.0.0 running GlassFish 4.1 Java 8 (Preconfigured - Docker)", "64bit Debian jessie v2.0.0 running GlassFish 4.0 Java 7 (Preconfigured - Docker)", "64bit Debian jessie v2.0.0 running Go 1.4 (Preconfigured - Docker)", "64bit Debian jessie v2.0.0 running Go 1.3 (Preconfigured - Docker)", "64bit Debian jessie v2.0.0 running Python 3.4 (Preconfigured - Docker)"])
end

retrieve_environment_info_1_patch = @patch function retrieve_environment_info(a...; b...)
    return Dict{String, Any}("EnvironmentInfo" => Any[Dict{String, Any}("SampleTimestamp" => "2015-08-20T22:23:17.703Z", "InfoType" => "tail", "Message" => "https://elasticbeanstalk-us-west-2-0123456789012.s3.amazonaws.com/resources/environments/logs/tail/e-fyqyju3yjs/i-09c1c867/TailLogs-1440109397703.out?AWSAccessKeyId=AKGPT4J56IAJ2EUBL5CQ&Expires=1440195891&Signature=n%2BEalOV6A2HIOx4Rcfb7LT16bBM%3D", "Ec2InstanceId" => "i-09c1c867")])
end

terminate_environment_1_patch = @patch function terminate_environment(a...; b...)
    return Dict{String, Any}("DateCreated" => "2015-08-12T18:52:53.622Z", "AbortableOperationInProgress" => false, "SolutionStackName" => "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 8 Java 8", "EndpointURL" => "awseb-e-f-AWSEBLoa-1I9XUMP4-8492WNUP202574.us-west-2.elb.amazonaws.com", "DateUpdated" => "2015-08-12T19:05:54.744Z", "Status" => "Terminating", "Tier" => Dict{String, Any}("Type" => "Standard", "Version" => " ", "Name" => "WebServer"), "EnvironmentName" => "my-env", "CNAME" => "my-env.elasticbeanstalk.com", "Health" => "Grey", "ApplicationName" => "my-app", "EnvironmentId" => "e-fh2eravpns")
end

update_application_1_patch = @patch function update_application(a...; b...)
    return Dict{String, Any}("Application" => Dict{String, Any}("Versions" => Any["2fba-stage-150819_234450", "bf07-stage-150820_214945", "93f8", "fd7c-stage-150820_000431", "22a0-stage-150819_185942"], "DateCreated" => "2015-08-13T19:15:50.449Z", "DateUpdated" => "2015-08-20T22:34:56.195Z", "Description" => "my Elastic Beanstalk application", "ApplicationName" => "my-app", "ConfigurationTemplates" => Any[]))
end

update_application_version_1_patch = @patch function update_application_version(a...; b...)
    return Dict{String, Any}("ApplicationVersion" => Dict{String, Any}("SourceBundle" => Dict{String, Any}("S3Key" => "my-app/22a0-stage-150819_185942.war", "S3Bucket" => "elasticbeanstalk-us-west-2-0123456789012"), "DateCreated" => "2015-08-19T18:59:17.646Z", "DateUpdated" => "2015-08-20T22:53:28.871Z", "Description" => "new description", "VersionLabel" => "22a0-stage-150819_185942", "ApplicationName" => "my-app"))
end

update_configuration_template_1_patch = @patch function update_configuration_template(a...; b...)
    return Dict{String, Any}("TemplateName" => "my-template", "DateCreated" => "2015-08-20T22:39:31Z", "DateUpdated" => "2015-08-20T22:43:11Z", "ApplicationName" => "my-app", "SolutionStackName" => "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 8 Java 8")
end

update_environment_1_patch = @patch function update_environment(a...; b...)
    return Dict{String, Any}("DateCreated" => "2015-02-03T23:04:54.453Z", "VersionLabel" => "v2", "SolutionStackName" => "64bit Amazon Linux running Tomcat 7", "EndpointURL" => "awseb-e-i-AWSEBLoa-1RDLX6TC9VUAO-0123456789.us-west-2.elb.amazonaws.com", "DateUpdated" => "2015-02-03T23:12:29.119Z", "Status" => "Updating", "Tier" => Dict{String, Any}("Type" => "Standard", "Version" => " ", "Name" => "WebServer"), "EnvironmentName" => "my-env", "CNAME" => "my-env.elasticbeanstalk.com", "Health" => "Grey", "ApplicationName" => "my-app", "EnvironmentId" => "e-szqipays4h")
end

update_environment_2_patch = @patch function update_environment(a...; b...)
    return Dict{String, Any}("DateCreated" => "2015-08-07T20:48:49.599Z", "AbortableOperationInProgress" => true, "VersionLabel" => "7f58-stage-150812_025409", "SolutionStackName" => "64bit Amazon Linux 2015.03 v2.0.0 running Tomcat 8 Java 8", "EndpointURL" => "awseb-e-w-AWSEBLoa-14XB83101Q4L-104QXY80921.sa-east-1.elb.amazonaws.com", "DateUpdated" => "2015-08-12T18:15:23.804Z", "Status" => "Updating", "Tier" => Dict{String, Any}("Type" => "Standard", "Version" => " ", "Name" => "WebServer"), "EnvironmentName" => "my-env", "CNAME" => "my-env.elasticbeanstalk.com", "Health" => "Grey", "ApplicationName" => "my-app", "EnvironmentId" => "e-wtp2rpqsej")
end

validate_configuration_settings_1_patch = @patch function validate_configuration_settings(a...; b...)
    return Dict{String, Any}("Messages" => Any[])
end


text:
# This test is derived from test cmcompute_naifid.c in Calceph version 2.3.2
# the test data files are copied from calceph-2.3.2.tar.gz
function testFunction2(testFile,testFile2,ephFiles,pflag)

    eph = Ephem(ephFiles)
    if pflag
        prefetch(eph)
    end

    con = constants(eph)
    AU = con[:AU]

    f = open(testFile);
    f2 = open(testFile2);
    for (ln,ln2) in zip(eachline(f),eachline(f2))
        elts=split(ln)
        elts2=split(ln2)
        jd0=parse(Float64,elts[1])
        @test parse(Float64,elts2[1]) == jd0
        target=parse(Int,elts2[2])
        center=parse(Int,elts2[3])
        targetold=parse(Int,elts[2])
        centerold=parse(Int,elts[3])
        dt = jd0 - trunc(Int,jd0)
        jd0 = trunc(Int,jd0) + 2.4515450000000000000E+06
        ref = [parse(Float64, x) for x in elts[4:end]]
        @test [parse(Float64, x) for x in elts2[4:end]] == ref

        if (target != naifId.id[:ttmtdb] && target != 15)
            for (unitold,ϵ) in [ (unitAU+unitDay,0.0),
                                 (unitKM+unitDay,0.0),
                                 (unitAU+unitSec,0.0),
                                 (unitKM+unitSec,0.0)]
               unit = unitold + useNaifId
               val = compute(eph, jd0, dt, target, center, unit)
               ref = compute(eph, jd0, dt, targetold, centerold, unitold)
               [(@test ref[i] ≈ val[i] atol=ϵ) for i in 1:6]
           end
       elseif (target == 15)
            targetN = 301
            for (unitold,ϵ) in [ (unitRad+unitDay,0.0),
                                 (unitRad+unitSec,0.0)]
               unit = unitold + useNaifId
               val = orient(eph, jd0, dt, targetN, unit)
               ref = compute(eph, jd0, dt, targetold, centerold, unitold)
               [(@test ref[i] ≈ val[i] atol=ϵ) for i in 1:6]
               val2 = orient(eph, jd0, dt, targetN, unit,3)
               @test length(val2) == 12
               [(@test val2[i] ≈ val[i] atol=ϵ) for i in 1:6]
           end
       elseif (target == naifId.id[:ttmtdb])
           for (unitold,ϵ) in [ (unitDay,0.0),
                                (unitSec,0.0)]
               unit = unitold + useNaifId
               val = compute(eph, jd0, dt, target, center, unit)
               ref = compute(eph, jd0, dt, targetold, centerold, unitold)
               @test ref[1] ≈ val[1] atol=ϵ
           end
       end

    end
    close(f)
    close(f2)
end

ext: jl


text:
# -*- coding: utf-8 -*-
# w32api_tests.jl

import W32API
using Base.Test

println("test W32API")
println("*** please look up MessageBox on your desktop ***")
println(@sprintf "[%08x]" W32API.MBox("-bodystring-", "titletest"))
println(@sprintf "[%08x]" W32API.MBox("bdy日本語漢字表示申能", "ttl表示申"))
println(@sprintf "[%08x]" W32API.MBoxA("-AbodyAstringA-", "AtitleAtestA"))
println(@sprintf "[%08x]" W32API.MBoxA("brokenOK漢表申能", "brokenOK表示申"))
wt = Array{UInt16, 1}([0x0074, 0x0074, 0x006c]) # w"ttl"
wb = Array{UInt16, 1}([0x002d, 0x0062, 0x0064, 0x0079, 0x002d]) # w"-bdy-"
println(@sprintf "[%08x]" W32API.MBoxW(wb, wt))
println(@sprintf "[%08x]" W32API.MBoxW("-callWbodyUTF8-", "callWtitleUTF8"))
wt = Array{UInt16, 1}([0x0074, 0x0074, 0x006c, 0x8868, 0x793a, 0x7533])
wb = Array{UInt16, 1}([0x0062, 0x0064, 0x0079,
  0x65e5, 0x672c, 0x8a9e, 0x6f22, 0x5b57, 0x8868, 0x793a, 0x7533, 0x80fd])
println(@sprintf "[%08x]" W32API.MBoxW(wb, wt))
println(@sprintf "[%08x]" W32API.MBoxW("callW表申能UTF8", "callW表示申UTF8"))
@test W32API.MBox("ok", "last") == 0x00000001
# @test W32API.MBox(wb, wt) != 0x00000001 # no function


text:
# This file is a part of Actors.jl, licensed under the MIT License (MIT).

import EasyPkg
EasyPkg.run_all_tests()


text: @show ARGS


text:
using Random
using CUDAdrv
using CUDAnative
using CuArrays
using StaticArrays
using GPUifyLoops # for @unroll macro

# Base.@propagate_inbounds function Base.getindex(v::MArray, i::Int)
#     @boundscheck checkbounds(v,i)
#     T = eltype(v)
#
#     @assert isbitstype(T) "Only isbits is supported on the device"
#     return GC.@preserve v begin
#         ptr  = pointer_from_objref(v)
#         dptr = reinterpret(CUDAnative.DevicePtr{T, CUDAnative.AS.Local}, ptr)
#         unsafe_load(dptr, i)
#     end
# end
#
# Base.@propagate_inbounds function Base.setindex!(v::MArray, val, i::Int)
#     @boundscheck checkbounds(v,i)
#     T = eltype(v)
#
#     @assert isbitstype(T) "Only isbits is supported on the device"
#     GC.@preserve v begin
#     ptr  = pointer_from_objref(v)
#         dptr = reinterpret(CUDAnative.DevicePtr{T, CUDAnative.AS.Local}, ptr)
#         unsafe_store!(dptr, convert(T, val), i)
#     end
#
#     return val
# end

for (f, T) in Base.Iterators.product((:add, :mul, :sub), (Float32, Float64))
    name = Symbol("$(f)_float_contract")
    if T === Float32
        llvmt = "float"
    elseif T === Float64
        llvmt = "double"
    end

    ir = """
    %x = f$f contract $llvmt %0, %1
    ret $llvmt %x
    """
    @eval begin
        # the @pure is necessary so that we can constant propagate.
        Base.@pure function $name(a::$T, b::$T)
            @Base._inline_meta
            Base.llvmcall($ir, $T, Tuple{$T, $T}, a, b)
        end
    end
end

# note the order of the fields below is also assumed in the code.
const _nstate = 5
const _ρ, _U, _V, _W, _E = 1:_nstate

const _nvgeo = 14
const _ξx, _ηx, _ζx, _ξy, _ηy, _ζy, _ξz, _ηz, _ζz, _MJ, _MJI,
       _x, _y, _z = 1:_nvgeo

Base.@irrational gdm1 0.4 BigFloat(0.4)

function volumerhs!(::Val{N}, rhs, Q, vgeo, D, nelem, t) where {N}
  Nq = N + 1

  s_D = @cuStaticSharedMem eltype(D) (Nq, Nq)
  s_F = @cuStaticSharedMem eltype(Q) (Nq, Nq, Nq, _nstate)
  s_G = @cuStaticSharedMem eltype(Q) (Nq, Nq, Nq, _nstate)
  s_H = @cuStaticSharedMem eltype(Q) (Nq, Nq, Nq, _nstate)

  e = blockIdx().x
  k = threadIdx().z
  j = threadIdx().y
  i = threadIdx().x

  @inbounds begin
    # fetch D into shared
    k == 1 && (s_D[i, j] = D[i, j])

    # Load values will need into registers
    MJ = vgeo[i, j, k, _MJ, e]
    ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
    ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
    ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]
    x, y, z = vgeo[i,j,k,_x,e], vgeo[i,j,k,_y,e], vgeo[i,j,k,_z,e]

    U, V, W = Q[i, j, k, _U, e], Q[i, j, k, _V, e], Q[i, j, k, _W, e]
    ρ, E = Q[i, j, k, _ρ, e], Q[i, j, k, _E, e]

    P = gdm1*(E - (U^2 + V^2 + W^2)/(2*ρ))

    ρinv = 1 / ρ

    fluxρ_x = U
    fluxU_x = ρinv * U * U + P
    fluxV_x = ρinv * U * V
    fluxW_x = ρinv * U * W
    fluxE_x = ρinv * U * (E + P)

    fluxρ_y = V
    fluxU_y = ρinv * V * U
    fluxV_y = ρinv * V * V + P
    fluxW_y = ρinv * V * W
    fluxE_y = ρinv * V * (E + P)

    fluxρ_z = W
    fluxU_z = ρinv * W * U
    fluxV_z = ρinv * W * V
    fluxW_z = ρinv * W * W + P
    fluxE_z = ρinv * W * (E + P)

    s_F[i, j, k, _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y + ξz * fluxρ_z)
    s_F[i, j, k, _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y + ξz * fluxU_z)
    s_F[i, j, k, _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y + ξz * fluxV_z)
    s_F[i, j, k, _W] = MJ * (ξx * fluxW_x + ξy * fluxW_y + ξz * fluxW_z)
    s_F[i, j, k, _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y + ξz * fluxE_z)

    s_G[i, j, k, _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y + ηz * fluxρ_z)
    s_G[i, j, k, _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y + ηz * fluxU_z)
    s_G[i, j, k, _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y + ηz * fluxV_z)
    s_G[i, j, k, _W] = MJ * (ηx * fluxW_x + ηy * fluxW_y + ηz * fluxW_z)
    s_G[i, j, k, _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y + ηz * fluxE_z)

    s_H[i, j, k, _ρ] = MJ * (ζx * fluxρ_x + ζy * fluxρ_y + ζz * fluxρ_z)
    s_H[i, j, k, _U] = MJ * (ζx * fluxU_x + ζy * fluxU_y + ζz * fluxU_z)
    s_H[i, j, k, _V] = MJ * (ζx * fluxV_x + ζy * fluxV_y + ζz * fluxV_z)
    s_H[i, j, k, _W] = MJ * (ζx * fluxW_x + ζy * fluxW_y + ζz * fluxW_z)
    s_H[i, j, k, _E] = MJ * (ζx * fluxE_x + ζy * fluxE_y + ζz * fluxE_z)

    r_rhsρ = pi*(-CUDAnative.sin(pi*t)*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 2*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y) +
                 CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 3*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*z) +
                 2*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 3*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 3*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z))
    r_rhsU = pi*(-10*CUDAnative.sin(pi*t)*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 15*CUDAnative.sin(pi*t)*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*z) -
                 15*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 10*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 5*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 30*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*z) -
                 3*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 15*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y) +
                 9*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 15*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 6*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y) +
                 18*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 2*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*z))*CUDAnative.cos(pi*y)/5
    r_rhsV = pi*(-10*CUDAnative.sin(pi*t)*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 15*CUDAnative.sin(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 3*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 9*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 10*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 5*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 6*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y) -
                 18*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 15*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y) +
                 15*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 15*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 2*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*z) +
                 30*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z))*CUDAnative.sin(pi*x)/5
    r_rhsW = pi*(-10*CUDAnative.sin(pi*t)*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 15*CUDAnative.sin(pi*t) -
                 15*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.sin(pi*(x + y))*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*z) +
                 15*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*z) +
                 36*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 18*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 4*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y) +
                 18*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 4*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y) +
                 30*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*z)*CUDAnative.cos(pi*(x + y)) -
                 2*CUDAnative.cos(pi*t))*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*y)/5
    r_rhsE = pi*(3*(1 - CUDAnative.cos(2*pi*x))*
                 (1 - CUDAnative.cos(4*pi*z))*(CUDAnative.cos(2*pi*t) + 1)*
                 (CUDAnative.cos(2*pi*y) + 1)/512 - 5*CUDAnative.sin(pi*t)*
                 CUDAnative.sin(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 4*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 8*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y) -
                 2*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 9*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 18*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 4*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 3*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 8*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 6*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 14*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 9*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 7*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*z)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y) -
                 18*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 7*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) -
                 700*CUDAnative.sin(pi*x)*CUDAnative.sin(pi*y)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*z) +
                 14*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 700*CUDAnative.sin(pi*x)*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z) +
                 700*CUDAnative.cos(pi*t)*CUDAnative.cos(pi*x)*CUDAnative.cos(pi*y)*CUDAnative.cos(pi*z))/5

    sync_threads()

    # loop of ξ-grid lines
    @unroll for n = 1:Nq
      Dni = s_D[n, i]
      Dnj = s_D[n, j]
      Dnk = s_D[n, k]

      r_rhsρ += Dni * s_F[n, j, k, _ρ]
      r_rhsρ += Dnj * s_G[i, n, k, _ρ]
      r_rhsρ += Dnk * s_H[i, j, n, _ρ]

      r_rhsU += Dni * s_F[n, j, k, _U]
      r_rhsU += Dnj * s_G[i, n, k, _U]
      r_rhsU += Dnk * s_H[i, j, n, _U]

      r_rhsV += Dni * s_F[n, j, k, _V]
      r_rhsV += Dnj * s_G[i, n, k, _V]
      r_rhsV += Dnk * s_H[i, j, n, _V]

      r_rhsW += Dni * s_F[n, j, k, _W]
      r_rhsW += Dnj * s_G[i, n, k, _W]
      r_rhsW += Dnk * s_H[i, j, n, _W]

      r_rhsE += Dni * s_F[n, j, k, _E]
      r_rhsE += Dnj * s_G[i, n, k, _E]
      r_rhsE += Dnk * s_H[i, j, n, _E]
    end

    MJI = vgeo[i, j, k, _MJI, e]

    rhs[i, j, k, _U, e] += MJI * r_rhsU
    rhs[i, j, k, _V, e] += MJI * r_rhsV
    rhs[i, j, k, _W, e] += MJI * r_rhsW
    rhs[i, j, k, _ρ, e] += MJI * r_rhsρ
    rhs[i, j, k, _E, e] += MJI * r_rhsE
  end

  nothing
end

function volumerhs_trigpi!(::Val{N}, rhs, Q, vgeo, D, nelem, t) where {N}
  Nq = N + 1

  s_D = @cuStaticSharedMem eltype(D) (Nq, Nq)
  s_F = @cuStaticSharedMem eltype(Q) (Nq, Nq, Nq, _nstate)
  s_G = @cuStaticSharedMem eltype(Q) (Nq, Nq, Nq, _nstate)
  s_H = @cuStaticSharedMem eltype(Q) (Nq, Nq, Nq, _nstate)

  e = blockIdx().x
  k = threadIdx().z
  j = threadIdx().y
  i = threadIdx().x

  @inbounds begin
    # fetch D into shared
    k == 1 && (s_D[i, j] = D[i, j])

    # Load values will need into registers
    MJ = vgeo[i, j, k, _MJ, e]
    ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
    ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
    ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]
    x, y, z = vgeo[i,j,k,_x,e], vgeo[i,j,k,_y,e], vgeo[i,j,k,_z,e]

    U, V, W = Q[i, j, k, _U, e], Q[i, j, k, _V, e], Q[i, j, k, _W, e]
    ρ, E = Q[i, j, k, _ρ, e], Q[i, j, k, _E, e]

    P = gdm1*(E - (U^2 + V^2 + W^2)/(2*ρ))

    ρinv = 1 / ρ

    fluxρ_x = U
    fluxU_x = ρinv * U * U + P
    fluxV_x = ρinv * U * V
    fluxW_x = ρinv * U * W
    fluxE_x = ρinv * U * (E + P)

    fluxρ_y = V
    fluxU_y = ρinv * V * U
    fluxV_y = ρinv * V * V + P
    fluxW_y = ρinv * V * W
    fluxE_y = ρinv * V * (E + P)

    fluxρ_z = W
    fluxU_z = ρinv * W * U
    fluxV_z = ρinv * W * V
    fluxW_z = ρinv * W * W + P
    fluxE_z = ρinv * W * (E + P)

    s_F[i, j, k, _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y + ξz * fluxρ_z)
    s_F[i, j, k, _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y + ξz * fluxU_z)
    s_F[i, j, k, _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y + ξz * fluxV_z)
    s_F[i, j, k, _W] = MJ * (ξx * fluxW_x + ξy * fluxW_y + ξz * fluxW_z)
    s_F[i, j, k, _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y + ξz * fluxE_z)

    s_G[i, j, k, _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y + ηz * fluxρ_z)
    s_G[i, j, k, _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y + ηz * fluxU_z)
    s_G[i, j, k, _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y + ηz * fluxV_z)
    s_G[i, j, k, _W] = MJ * (ηx * fluxW_x + ηy * fluxW_y + ηz * fluxW_z)
    s_G[i, j, k, _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y + ηz * fluxE_z)

    s_H[i, j, k, _ρ] = MJ * (ζx * fluxρ_x + ζy * fluxρ_y + ζz * fluxρ_z)
    s_H[i, j, k, _U] = MJ * (ζx * fluxU_x + ζy * fluxU_y + ζz * fluxU_z)
    s_H[i, j, k, _V] = MJ * (ζx * fluxV_x + ζy * fluxV_y + ζz * fluxV_z)
    s_H[i, j, k, _W] = MJ * (ζx * fluxW_x + ζy * fluxW_y + ζz * fluxW_z)
    s_H[i, j, k, _E] = MJ * (ζx * fluxE_x + ζy * fluxE_y + ζz * fluxE_z)

    r_rhsρ = pi*(-CUDAnative.sinpi(t)*CUDAnative.sinpi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 2*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y) +
                 CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 3*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(z) +
                 2*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 3*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 3*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z))
    r_rhsU = pi*(-10*CUDAnative.sinpi(t)*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 15*CUDAnative.sinpi(t)*CUDAnative.sinpi(x)*CUDAnative.cospi(z) -
                 15*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 10*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 5*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 30*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(z) -
                 3*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 15*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y) +
                 9*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 15*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 6*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y) +
                 18*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 2*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(z))*CUDAnative.cospi(y)/5
    r_rhsV = pi*(-10*CUDAnative.sinpi(t)*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 15*CUDAnative.sinpi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 3*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 9*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 10*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 5*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 6*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y) -
                 18*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 15*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y) +
                 15*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 15*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 2*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(z) +
                 30*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z))*CUDAnative.sinpi(x)/5
    r_rhsW = pi*(-10*CUDAnative.sinpi(t)*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 15*CUDAnative.sinpi(t) -
                 15*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.sinpi((x + y))*CUDAnative.cospi(t)*CUDAnative.cospi(z) +
                 15*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(z) +
                 36*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 18*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 4*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y) +
                 18*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 4*CUDAnative.cospi(t)*CUDAnative.cospi(y) +
                 30*CUDAnative.cospi(t)*CUDAnative.cospi(z)*CUDAnative.cospi((x + y)) -
                 2*CUDAnative.cospi(t))*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(y)/5
    r_rhsE = pi*(3*(1 - CUDAnative.cos(2*pi*x))*
                 (1 - CUDAnative.cos(4*pi*z))*(CUDAnative.cos(2*pi*t) + 1)*
                 (CUDAnative.cos(2*pi*y) + 1)/512 - 5*CUDAnative.sinpi(t)*
                 CUDAnative.sinpi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 4*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 8*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y) -
                 2*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 9*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 18*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 4*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 3*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 8*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 6*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 14*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 9*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 7*CUDAnative.sinpi(x)*CUDAnative.sinpi(z)*CUDAnative.cospi(t)*CUDAnative.cospi(y) -
                 18*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 7*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) -
                 700*CUDAnative.sinpi(x)*CUDAnative.sinpi(y)*CUDAnative.cospi(t)*CUDAnative.cospi(z) +
                 14*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 700*CUDAnative.sinpi(x)*CUDAnative.cospi(t)*CUDAnative.cospi(y)*CUDAnative.cospi(z) +
                 700*CUDAnative.cospi(t)*CUDAnative.cospi(x)*CUDAnative.cospi(y)*CUDAnative.cospi(z))/5

    sync_threads()

    # loop of ξ-grid lines
    @unroll for n = 1:Nq
      Dni = s_D[n, i]
      Dnj = s_D[n, j]
      Dnk = s_D[n, k]

      r_rhsρ += Dni * s_F[n, j, k, _ρ]
      r_rhsρ += Dnj * s_G[i, n, k, _ρ]
      r_rhsρ += Dnk * s_H[i, j, n, _ρ]

      r_rhsU += Dni * s_F[n, j, k, _U]
      r_rhsU += Dnj * s_G[i, n, k, _U]
      r_rhsU += Dnk * s_H[i, j, n, _U]

      r_rhsV += Dni * s_F[n, j, k, _V]
      r_rhsV += Dnj * s_G[i, n, k, _V]
      r_rhsV += Dnk * s_H[i, j, n, _V]

      r_rhsW += Dni * s_F[n, j, k, _W]
      r_rhsW += Dnj * s_G[i, n, k, _W]
      r_rhsW += Dnk * s_H[i, j, n, _W]

      r_rhsE += Dni * s_F[n, j, k, _E]
      r_rhsE += Dnj * s_G[i, n, k, _E]
      r_rhsE += Dnk * s_H[i, j, n, _E]
    end

    MJI = vgeo[i, j, k, _MJI, e]

    rhs[i, j, k, _U, e] += MJI * r_rhsU
    rhs[i, j, k, _V, e] += MJI * r_rhsV
    rhs[i, j, k, _W, e] += MJI * r_rhsW
    rhs[i, j, k, _ρ, e] += MJI * r_rhsρ
    rhs[i, j, k, _E, e] += MJI * r_rhsE
  end

  nothing
end

function main()
    N = 4
    nelem = 4000
    ntrials = 1
    DFloat = Float32

    rnd = MersenneTwister(0)

    Nq = N + 1
    Q = 1 .+ CuArray(rand(rnd, DFloat, Nq, Nq, Nq, _nstate, nelem))
    Q[:, :, :, _E, :] .+= 20
    vgeo = CuArray(rand(rnd, DFloat, Nq, Nq, Nq, _nvgeo, nelem))

    # Make sure the entries of the mass matrix satisfy the inverse relation
    vgeo[:, :, :, _MJ, :] .+= 3
    vgeo[:, :, :, _MJI, :] .= 1 ./ vgeo[:, :, :, _MJ, :]

    D = CuArray(rand(rnd, DFloat, Nq, Nq))

    rhs = CuArray(zeros(DFloat, Nq, Nq, Nq, _nstate, nelem))

    t = DFloat(1)
    @cuda(threads=(N+1, N+1, N+1), blocks=nelem, maxregs=255,
          volumerhs!(Val(N), rhs, Q, vgeo, D, nelem, t))

    CUDAdrv.@profile  begin
        @cuda(threads=(N+1, N+1, N+1), blocks=nelem, maxregs=255,
              volumerhs!(Val(N), rhs, Q, vgeo, D, nelem, t))
    end

    @cuda(threads=(N+1, N+1, N+1), blocks=nelem, maxregs=255,
          volumerhs_trigpi!(Val(N), rhs, Q, vgeo, D, nelem, t))

    CUDAdrv.@profile  begin
        @cuda(threads=(N+1, N+1, N+1), blocks=nelem, maxregs=255,
              volumerhs_trigpi!(Val(N), rhs, Q, vgeo, D, nelem, t))
    end

    #=
    CUDAnative.@device_code dir="jlir" begin
        @cuda(threads=(N+1, N+1, N+1), blocks=nelem, maxregs=255,
              volumerhs!(Val(N), rhs, Q, vgeo, D, nelem, t))
    end
    =#

    nothing
end

main()


text:
#Script to make a series of demonstration plots
#27/07/2017

using Plots, JLD, LaTeXStrings
pyplot()
include("functions.jl")
PyPlot.close("all")

#Input information
atom = "87Rb"
nn = 50
ll = 0
jj = 0.5
mj = 0.5

alpha_c = getalpha(atom)

#Calculate wave function
normY_sol, rr = numerovfunc(atom,nn,ll,jj)

#Rescale for plotting
plotscale = sqrt(rr)
probamp = (normY_sol.*plotscale).^2

PyPlot.figure()
if nn>20
  plot(plotscale[plotscale .> sqrt(alpha_c^(1/3))],normY_sol[plotscale .> sqrt(alpha_c^(1/3))],linewidth=2)
else
  plot(plotscale,normY_sol,linewidth=2)
end
plot!(xlabel=L"(r/a_0)^{1/2}")
plot!(ylabel=L"r^{1/2}R(r) \,(a_0^{-1})")
plot!(title=string(atom," radial wavefunction |n,l,j⟩ = |",string(nn),",",string(ll),",",string(jj),"⟩"))
plot!(leg=false)
gui()

PyPlot.figure()
if nn>20
  plot(rr[rr .> alpha_c^(1/3)],probamp[rr .> alpha_c^(1/3)],linewidth=2)
else
  plot(rr,probamp,linewidth=2)
end
plot!(xlabel=L"r/a_0")
plot!(ylabel=L"|rR(r)|^2 \,(a_0^{-1})")
plot!(title=string(atom," radial probability density |n,l,j⟩ = |",string(nn),",",string(ll),",",string(jj),"⟩"))
plot!(leg=false)
gui()

#Compare groundstate to self-consistent field method
scfdict = load("Data/SCFcalc.jld")
plotscale2 = scfdict["plotscale2"]
sqrtrR = scfdict["sqrtrR"]
normY_solg, rrg = numerovfunc(atom,5,0,0.5)

PyPlot.figure()
plot(sqrt(rrg),-normY_solg,linewidth=2,label="Numerov")
scatter!(plotscale2,sqrtrR,label="Callaway")
plot!(xlabel=L"r/a_0")
plot!(ylabel=L"|r^{1/2}R(r)|^2 \,(a_0^{-1})")
plot!(title="Rubidium Groundstate Radial Wavefunction")
gui()

#Lifetimes
ltdict = load("Data/Lifetimes.jld")
sexpt = ltdict["sexpt"]
sexpterr = ltdict["sexpterr"]
pexpt = ltdict["pexpt"]
pexpterr = ltdict["pexpterr"]
dexpt = ltdict["dexpt"]
dexpterr = ltdict["dexpterr"]

#S1_2
nns = collect(28:45)
calcs = zeros(length(nns))
for kk = 1:length(nns)
  calcs[kk] = Radiative_lifetimes(atom,nns[kk],0,0.5)
end
calcs = calcs*1e6

#P3_2
nnp = collect(34:44)
calcp = zeros(length(nnp))
for kk = 1:length(nnp)
  calcp[kk] = Radiative_lifetimes(atom,nnp[kk],1,1.5)
end
calcp = calcp*1e6

#D5_2
nnd = collect(29:44)
calcd = zeros(length(nnd))
for kk = 1:length(nnd)
  calcd[kk] = Radiative_lifetimes(atom,nnd[kk],2,2.5)
end
calcd = calcd*1e6

PyPlot.figure()
scatter(nns,calcs,label="Calculated")
scatter!(nns,sexpt[:],yerr=repmat(sexpterr[:],1,2),label="Experimental")
plot!(xlabel=L"n")
plot!(ylabel=L"τ\, (μ\mathrm{s})")
plot!(title="Plot of calculated and experimental lifetimes for \n |n,l,j⟩=|n,0,0.5⟩ states")
gui()

PyPlot.figure()
scatter(nnp,calcp,label="Calculated")
scatter!(nnp,pexpt[:],yerr=repmat(pexpterr[:],1,2),label="Experimental")
plot!(xlabel=L"n")
plot!(ylabel=L"τ\, (μ\mathrm{s})")
plot!(title="Plot of calculated and experimental lifetimes for \n |n,l,j⟩=|n,1,1.5⟩ states")
gui()

PyPlot.figure()
scatter(nnd,calcd,label="Calculated")
scatter!(nnd,dexpt[:],yerr=repmat(dexpterr[:],1,2),label="Experimental")
plot!(xlabel=L"n")
plot!(ylabel=L"τ\, (μ\mathrm{s})")
plot!(title="Plot of calculated and experimental lifetimes for \n |n,l,j⟩=|n,2,2.5⟩ states")
gui()

#Blockade shift
RRSI, θ, blockadeshiftmeshGHz, C_6val = BlockadeShift(atom,nn,ll,jj,mj)

PyPlot.figure()
heatmap(RRSI*1e6,θ,blockadeshiftmeshGHz*1e3)
plot!(xlabel=L"R\, (μ\mathrm{m})")
plot!(ylabel=L"θ")
plot!(zlabel=L"ΔW\, (\mathrm{MHz})")
plot!(title=string("Calculated Rydberg blockade shift for \n |n,l,j,mⱼ⟩ = |",string(nn),",",string(ll),",",string(jj),",",string(mj),"⟩ state of ",atom," [MHz]"))
gui()

#C₆
c6dict = load("Data/C6data.jld")
nnc6 = c6dict["nnc6"]
#Convert from GHz⋅m⁻⁶ to GHz⋅μm⁻⁶
c6S1_2 = c6dict["c6S1_2"]*1e36
c6D3_2 = c6dict["c6D3_2"]*1e36
c6D5_2 = c6dict["c6D5_2"]*1e36

PyPlot.figure()
plot(nnc6,abs(c6S1_2),linewidth=2,marker=(:circle))
plot!(yaxis=(:log))
plot!(xlabel=L"n")
plot!(ylabel=L"|C₆|\, (\mathrm{GHz}⋅μ\mathrm{m}⁻⁶)")
plot!(title="Calculated C₆ for |n,l,j⟩ = |n,0,0.5⟩")
plot!(leg=false)
gui()

PyPlot.figure()
plot(nnc6,abs(c6D3_2),linewidth=2,marker=(:circle))
plot!(yaxis=(:log))
plot!(xlabel=L"n")
plot!(ylabel=L"|C₆|\, (\mathrm{GHz}⋅μ\mathrm{m}⁻⁶)")
plot!(title="Calculated C₆ for |n,l,j⟩ = |n,2,1.5⟩")
plot!(leg=false)
gui()

PyPlot.figure()
plot(nnc6,abs(c6D5_2),linewidth=2,marker=(:circle))
plot!(yaxis=(:log))
plot!(xlabel=L"n")
plot!(ylabel=L"|C₆|\, (\mathrm{GHz}⋅μ\mathrm{m}⁻⁶)")
plot!(title="Calculated C₆ for |n,l,j⟩ = |n,2,2.5⟩")
plot!(leg=false)
gui()

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
include("tcm_estimateForest.jl")

function tcm_estimateForest_leafLimit(N, assortments, transaction_counts, time_limit_overall, perTol, leaf_limit, initial_left, initial_right, initial_product, initial_isLeaf)
    M = size(assortments, 1)

    time_limit_EM = 60.0

    start_time = time()

    forest_left = deepcopy(initial_left)
    forest_right = deepcopy(initial_right)
    forest_product = deepcopy(initial_product)
    forest_isLeaf = deepcopy(initial_isLeaf)

    K = length(forest_left)

    #
    initial_lambda = 1/K * ones(K)
    lambda, v, loglik, elapsed_time = tcm_forest_EM(N, K, forest_left, forest_right, forest_product, forest_isLeaf, assortments, transaction_counts, time_limit_EM, initial_lambda)

    # @show loglik

    u = transaction_counts ./ v
    u[ isinf.(u)] = 0.0
    u[ isnan.(u)] = 0.0

    # NB - re. isinf.(...), this should not be necessary, if the initial model is such that every item assortment pair has a positive choice probability;
    # this is guaranteed for the ranking model because it starts from a set of rankings corresponding to an independent demand model. If the forest
    # is initialized with the ranking model, then the same condition will hold.
    # for isnan.(...), this is definitely necessary, because transaction_counts and v will both be zero at items that are not part of the assortment.

    # perTol = 1e-7
    iter = 0;

    println("CG iter ", iter, " -- time elapsed = ", 0.0, " log likelihood = ",loglik)

    while (time() - start_time < time_limit_overall)
        iter += 1

        u = transaction_counts ./ v
        u[ isinf.(u)] = 0.0
        u[ isnan.(u)] = 0.0
        single_left, single_right, single_product, single_isLeaf, objval = tcm_solveTreeSubproblem_leafLimit(N, u, assortments, leaf_limit)
        # error("Stop here")

        #orderings = [orderings; single_ordering']
        push!(forest_left, single_left)
        push!(forest_right, single_right)
        push!(forest_product, single_product)
        push!(forest_isLeaf, single_isLeaf)

        K += 1

        initial_fudge = 0.999
        initial_lambda = [ initial_fudge*lambda; (1-initial_fudge)]
        lambda, v, new_loglik, elapsed_time = tcm_forest_EM(N, K, forest_left, forest_right, forest_product, forest_isLeaf, assortments, transaction_counts, time_limit_EM, initial_lambda)

        println("CG iter ", iter, " -- time elapsed = ", time() - start_time, " log likelihood = ",new_loglik)

        if ( (new_loglik - loglik) / abs(loglik) < perTol )
            temp = (new_loglik - loglik) / abs(loglik)
            @show temp
            @show new_loglik
            @show loglik

            loglik = new_loglik
            break
        end

        loglik = new_loglik

    end

    end_time = time()

    elapsed_time = end_time - start_time

    return lambda, forest_left, forest_right, forest_product, forest_isLeaf, loglik, elapsed_time
end

function tcm_solveTreeSubproblem_leafLimit(N, u, assortments, leaf_limit)
    M = size(assortments,1)

    single_left = Int64[0]
    single_right = Int64[0]
    single_product = Int64[N]
    single_isLeaf = Bool[true]
    single_depth = Int64[1]
    single_include_set = Array{Int64}[ [N] ];
    single_assortment_binary = Array{Int64}[ ones(M) ]
    single_reduced_cost = Float64[0.0];
    single_isClosed = Bool[false]

    iter = 0 ;

    while (true)
        current_reduced_cost = sum( single_reduced_cost[single_isLeaf])
        println("\t\t Subproblem Iter : ", iter, " -- reduced cost: ", current_reduced_cost, " -- num leaves: ", sum(single_isLeaf))
        candidate_leaves = find(single_isLeaf .&  (.!single_isClosed) )

        if (isempty(candidate_leaves))
            println("\t\t\t No more eligible leaves to split; exiting procedure...")
            break;
        end

        if (sum(single_isLeaf) >= leaf_limit)
            println("\t\t\t Leaf limit reached; exiting procedure...")
            break;
        end

        num_candL = length(candidate_leaves)
        candidate_reduced_costs = -Inf* ones(num_candL)
        best_split_p_by_leaf = zeros(Int64, num_candL)
        best_left_leaf_p_by_leaf  = zeros(Int64, num_candL)
        best_right_leaf_p_by_leaf  = zeros(Int64, num_candL)
        best_left_reduced_cost_by_leaf  = zeros(Float64, num_candL)
        best_right_reduced_cost_by_leaf  = zeros(Float64, num_candL)
        best_left_assortment_binary_by_leaf = Array{Array{Int64}}(num_candL)
        best_right_assortment_binary_by_leaf = Array{Array{Int64}}(num_candL)

        for ell_ind in 1:num_candL

            leaf = candidate_leaves[ell_ind]

            assortment_binary = copy(single_assortment_binary[leaf])

            candidate_split_products = setdiff(collect(1:(N-1)), single_include_set[leaf])

            for split_p in candidate_split_products

                candidate_left_leaf_products = union(single_include_set[leaf], split_p)
                best_left_reduced_cost = -Inf
                best_left_p = 0;
                best_left_assortment_binary = zeros(M)
                for left_p in candidate_left_leaf_products
                    left_assortment_binary = assortment_binary .* assortments[:, split_p][:]
                    temp = dot( left_assortment_binary, u[:, left_p])
                    if (temp > best_left_reduced_cost)
                        best_left_reduced_cost = temp
                        best_left_p = left_p
                        best_left_assortment_binary = copy(left_assortment_binary)
                    end
                end

                candidate_right_leaf_products = copy(single_include_set[leaf])
                best_right_reduced_cost = -Inf
                best_right_p = 0;
                best_right_assortment_binary = zeros(M)
                for right_p in candidate_right_leaf_products
                    right_assortment_binary = assortment_binary .* (1 - assortments[:,split_p][:])
                    temp = dot(right_assortment_binary, u[:, right_p])
                    if (temp > best_right_reduced_cost)
                        best_right_reduced_cost = temp
                        best_right_p = right_p
                        best_right_assortment_binary = copy(right_assortment_binary)
                    end
                end

                new_reduced_cost = best_left_reduced_cost + best_right_reduced_cost

                if (new_reduced_cost > candidate_reduced_costs[ell_ind])
                    candidate_reduced_costs[ell_ind] = new_reduced_cost
                    best_split_p_by_leaf[ell_ind] = split_p
                    best_left_leaf_p_by_leaf[ell_ind] = best_left_p
                    best_right_leaf_p_by_leaf[ell_ind] = best_right_p
                    best_left_reduced_cost_by_leaf[ell_ind] = best_left_reduced_cost
                    best_right_reduced_cost_by_leaf[ell_ind] = best_right_reduced_cost
                    best_left_assortment_binary_by_leaf[ell_ind] = best_left_assortment_binary # VVM: should there be a copy here?
                    best_right_assortment_binary_by_leaf[ell_ind] = best_right_assortment_binary # VVM: should there be a copy here?
                end
            end
        end

        incumbent_reduced_costs = single_reduced_cost[candidate_leaves]

        # Find most improving leaf:
        ell_ind = indmax(candidate_reduced_costs - incumbent_reduced_costs)
        tol_RC = 1e-4

        if (incumbent_reduced_costs[ell_ind] + tol_RC < candidate_reduced_costs[ell_ind])
            # Splitting this leaf does improve things, so go ahead with it.

            leaf = candidate_leaves[ell_ind]

            # Grow tree.
            numNodes = length(single_left)
            push!(single_left, 0, 0)
            push!(single_right, 0, 0)
            push!(single_isLeaf, true, true)
            push!(single_product, best_left_leaf_p_by_leaf[ell_ind], best_right_leaf_p_by_leaf[ell_ind])
            push!(single_depth, single_depth[leaf]+1, single_depth[leaf]+1)
            push!(single_include_set, union(single_include_set[leaf],best_split_p_by_leaf[ell_ind]), copy(single_include_set[leaf]) )
            push!(single_assortment_binary, copy(best_left_assortment_binary_by_leaf[ell_ind]), copy(best_right_assortment_binary_by_leaf[ell_ind]))
            push!(single_reduced_cost, best_left_reduced_cost_by_leaf[ell_ind], best_right_reduced_cost_by_leaf[ell_ind])
            # push!(single_isClosed, false, false)

            # Update the leaf
            single_left[leaf] = numNodes+1
            single_right[leaf] = numNodes+2
            single_product[leaf] = best_split_p_by_leaf[ell_ind]
            # single_isClosed[leaf] = true
            single_isLeaf[leaf] = false
        else
            # Leaf does not improve
            println("\t\t\t No more improvement; exiting procedure...")
            break;
        end
    end

    optimal_reduced_cost = sum( single_reduced_cost[single_isLeaf] )

    return single_left, single_right, single_product, single_isLeaf, optimal_reduced_cost

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# returns the variable (or matrix), lagged p times,
# with the first p rows filled with ones (to avoid divide errors)
# remember to drop those rows before doing analysis
function lag(x,p)
	n = size(x,1)
    k = size(x,2)
	lagged_x = [ones(p,k); x[1:n-p,:]]
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
module PythonOT

using PyCall: PyCall

export emd,
    emd2,
    emd_1d,
    emd2_1d,
    sinkhorn,
    sinkhorn2,
    barycenter,
    barycenter_unbalanced,
    sinkhorn_unbalanced,
    sinkhorn_unbalanced2

const pot = PyCall.PyNULL()

include("lib.jl")
include("smooth.jl")

function __init__()
    return copy!(pot, PyCall.pyimport_conda("ot", "pot", "conda-forge"))
end

end


text:
# # HELPERS FOR `FIT` AND `UPDATE`

# either retrains new machine on all data or returns the machine
# trained using controls (embedded in `ic_model`):
function production_mach(iterated_model,
                         ic_model,
                         verbosity,
                         data...)
    train_mach = IterationControl.expose(ic_model)
    if iterated_model.retrain &&
        iterated_model.resampling !== nothing
        iteration_param = ic_model.iteration_parameter
        clone = deepcopy(iterated_model.model)
        # set iteration parameter to value at end of controlled
        # training:
        n_iterations = rget(train_mach.model, iteration_param)
        rset!(clone, iteration_param, n_iterations)
        prod_mach = machine(clone,
                            data...;
                            cache=iterated_model.cache)
        verbosity < 1 ||
            @info "Retraining on all provided data. "*
            "To suppress, specify `retrain=false`. "
        fit!(prod_mach, verbosity=verbosity - 1)
    else
        train_mach
    end
end

function finish(prod_mach, ic_model, iterated_model, creport, verbosity)

    iteration_param = ic_model.iteration_parameter

    verbosity < 1 ||
        @info "Total of $(rget(prod_mach.model, iteration_param)) "*
        "iterations. "

    fitresult = prod_mach
    cache = (ic_model=ic_model,
             iterated_model=deepcopy(iterated_model))

    report = (model_report=MLJBase.report(prod_mach),
              controls=creport,
              n_iterations=rget(prod_mach, iteration_param))
    return fitresult, cache, report
end

# # IMPLEMENTATION OF MLJ MODEL INTERFACE

function MLJBase.fit(iterated_model::EitherIteratedModel, verbosity, data...)

    model = deepcopy(iterated_model.model)

     # get name of iteration parameter:
     _iter = MLJBase.iteration_parameter(model)
    iteration_param = _iter === nothing ?
        iterated_model.iteration_parameter : _iter

    # instantiate `train_mach`:
    mach = if iterated_model.resampling === nothing
        machine(model, data...; cache=iterated_model.cache)
    else
        resampler = MLJBase.Resampler(model=model,
                                  resampling=iterated_model.resampling,
                                  measure=iterated_model.measure,
                                  weights=iterated_model.weights,
                                  class_weights=iterated_model.class_weights,
                                  operation=iterated_model.operation,
                                  check_measure=iterated_model.check_measure,
                                  cache=iterated_model.cache)
        machine(resampler, data..., cache=false)
    end

    # instantiate the object to be iterated using IterationControl.jl:
    ic_model = ICModel(mach, iteration_param, verbosity)

    # train with controls:
    creport = IterationControl.train!(ic_model,
                                      iterated_model.controls...,
                                      verbosity=verbosity)

    # retrains on all data if necessary:
    prod_mach = production_mach(iterated_model,
                                ic_model,
                                verbosity,
                                data...)

    return finish(prod_mach, ic_model, iterated_model, creport, verbosity)

end

function MLJBase.update(iterated_model::EitherIteratedModel,
                        verbosity,
                        old_fitresult,
                        old_cache,
                        data...)

    ic_model, old_iterated_model = old_cache

    # cold restart if anything but `controls` or `model` altered:
    if !MLJBase.is_same_except(iterated_model,
                              old_iterated_model,
                              :model, :controls)
        return fit(iterated_model, verbosity, data...)
    end

    # otherwise, continue training with existing `ic_model`:
    # train with controls:
    creport = IterationControl.train!(ic_model,
                                      iterated_model.controls...,
                                      verbosity=verbosity)

    # retrains on all data if necessary:
    prod_mach = production_mach(iterated_model,
                                ic_model,
                                verbosity,
                                data...)

    return finish(prod_mach, ic_model, iterated_model, creport, verbosity)

end

MLJBase.fitted_params(::EitherIteratedModel, fitresult) = (machine=fitresult,)

MLJBase.predict(::EitherIteratedModel, fitresult, Xnew) =
    predict(fitresult, Xnew)

MLJBase.transform(::EitherIteratedModel, fitresult, Xnew) =
    transform(fitresult, Xnew)

max_stars_repo_name: ablaom/MLJIteration

max_issues_repo_name: ablaom/MLJIteration

max_forks_repo_name: ablaom/MLJIteration


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.test.behaviour.graql.language.match;
#
# import typedb.core.test.behaviour.BehaviourTest;
# import io.cucumber.junit.Cucumber;
# import io.cucumber.junit.CucumberOptions;
# import org.junit.runner.RunWith;
#
# @RunWith(Cucumber.class)
# @CucumberOptions(
#         strict = true,
#         plugin = "pretty",
#         glue = "typedb.client.test.behaviour",
#         features = "external/typedblabs_behaviour/graql/language/match.feature",
#         tags = "not @ignore and not @ignore-client-java"
# )
# public class MatchTest extends BehaviourTest {
#     // ATTENTION:
#     // When you click RUN from within this class through Intellij IDE, it will fail.
#     // You can fix it by doing:
#     //
#     // 1) Go to 'Run'
#     // 2) Select 'Edit Configurations...'
#     // 3) Select 'Bazel test MatchTest'
#     //
#     // 4) Ensure 'Target Expression' is set correctly:
#     //    a) Use '//<this>/<package>/<name>:test-core' to test against typedb-core
#     //    b) Use '//<this>/<package>/<name>:test-kgms' to test against typedb-kgms
#     //
#     // 5) Update 'Bazel Flags':
#     //    a) Remove the line that says: '--test_filter=typedb.client.*'
#     //    b) Use the following Bazel flags:
#     //       --cache_test_results=no : to make sure you're not using cache
#     //       --test_output=streamed : to make sure all output is printed
#     //       --subcommands : to print the low-level commands and execution paths
#     //       --sandbox_debug : to keep the sandbox not deleted after test runs
#     //       --spawn_strategy=standalone : if you're on Mac, tests need permission to access filesystem (to run TypeDB)
#     //
#     // 6) Hit the RUN button by selecting the test from the dropdown menu on the top bar    private static TypeDBRunner runner;
# }


text:
cd(abs(parent( Path(string(@__FILE__)) ))) do
    @testset "Simple Path Usage" begin
        reg = "../src/FilePaths.jl"
        @test ispath(reg)

        p = Path(reg)

        @test p == p"../src/FilePaths.jl"
        @test string(p) == reg
        @test string(cwd()) == pwd()
        @test string(home()) == homedir()

        @test parts(p) == ("..", "src", "FilePaths.jl")
        @test hasparent(p)
        @test parent(p) == p"../src"
        @test parents(p) == [p"..", p"../src"]
        @test_throws ErrorException parents(p".")

        @test basename(p) == "FilePaths.jl"
        @test joinpath(parent(p), Path(basename(p))) == p
        @test filename(p) == "FilePaths"

        @test extension(p) == "jl"
        @test extension(p"../REQUIRE") == ""
        @test extensions(p"foo.tar.gz") == ["tar", "gz"]
        @test length(extensions(p"../REQUIRE")) == 0

        @test exists(p)
        @test !isabs(p)
        @test string(norm(p"../src/../src/FilePaths.jl")) == normpath("../src/../src/FilePaths.jl")
        @test string(abs(p)) == abspath(string(p))
        @test string(relative(p, home())) == relpath(string(p), homedir())
        @test uri(PosixPath("/foo/bar")) == URI("file:///foo/bar")
        @test_throws ErrorException uri(p"foo/bar")

        @test p"../src/FilePaths.jl" in glob(p"../src", "*.jl")

        s = stat(p)
        lstat(p)

        io = IOBuffer()
        show(io, s)
        show_str = copy(takebuf_string(io))
        #@test "device" in show_str
        #@test "blocks" in show_str

        @test size(p) == stat(p).size
        @test modified(p) == stat(p).mtime
        @test created(p) == stat(p).ctime

        @test isfile(p)
        @test isdir(parent(p))
        @test !islink(p)
        @test !issocket(p)
        @test !isfifo(p)
        @test !ischardev(p)
        @test !isblockdev(p)
    end
end

mktmpdir() do d
    cd(d) do
        @testset "Modifying Path Usage" begin
            new_path = p"foo"

            mkdir(new_path)
            @test_throws ErrorException mkdir(new_path)
            remove(new_path)

            new_path = p"foo/bar"
            @test_throws ErrorException mkdir(new_path)
            mkdir(new_path; recursive=true)
            mkdir(new_path; recursive=true, exist_ok=true)

            other_path = p"car/bar"
            copy(new_path, other_path; recursive=true)
            copy(new_path, other_path; exist_ok=true, overwrite=true)
            @test_throws ErrorException copy(new_path, other_path)
            @test_throws ErrorException copy(p"badpath", other_path; exist_ok=true, overwrite=true)
            remove(p"car"; recursive=true)

            move(new_path, other_path; recursive=true)
            mkdir(new_path; recursive=true)
            move(new_path, other_path; exist_ok=true, overwrite=true)
            @test_throws ErrorException move(new_path, other_path)
            @test_throws ErrorException move(p"badpath", other_path; exist_ok=true, overwrite=true)
            remove(p"car"; recursive=true)

            mkdir(new_path; recursive=true)

            symlink(new_path, p"mysymlink")
            symlink(new_path, p"mysymlink"; exist_ok=true, overwrite=true)
            @test_throws ErrorException symlink(new_path, p"mysymlink")
            @test_throws ErrorException symlink(p"badpath", p"mysymlink"; exist_ok=true, overwrite=true)

            touch(p"newfile")
            mktmp(d) do f, io
                println(f)
                println(io)
            end

            @unix_only if ENV["USER"] =="root"
                chown(p"newfile", "nobody", "nogroup"; recursive=true)
            else
                @test_throws ErrorException chown(p"newfile", "nobody", "nogroup"; recursive=true)
            end

            @windows_only @test_throws ErrorException chown(p"newfile", "nobody", "nogroup"; recursive=true)

            chmod(p"newfile", user=(READ+WRITE+EXEC), group=(READ+EXEC), other=READ)
            @test string(mode(p"newfile")) == "-rwxr-xr--"
            @test isexecutable(p"newfile")
            @test iswritable(p"newfile")
            @test isreadable(p"newfile")

            chmod(p"newfile", "-x")
            @test !isexecutable(p"newfile")

            @test string(mode(p"newfile")) == "-rw-r--r--"
            chmod(p"newfile", "+x")
            write(p"newfile", "foobar")
            @test read(p"newfile") == "foobar"
            chmod(p"newfile", "u=rwx")

            chmod(new_path, mode(p"newfile"); recursive=true)

        end
    end
end


text:
export MDFFile, MDFFileV1, MDFFileV2, addTrailingSingleton, addLeadingSingleton

abstract type MDFFile <: MPIFile end

# We use a dedicated type for v1 and v2. If both versions
# are the same we use the abstract type MDFFile
mutable struct MDFFileV1{T} <: MDFFile
  filename::String
  file::HDF5.File
  mmap_measData::T
end

function MDFFileV1(filename::String, file=h5open(filename,"r"))
  tdExists = haskey(file, "/measurement/dataTD")
  if haskey(file, "/measurement/dataTD")
    mmap_measData = HDF5.readmmap(file["/measurement/dataTD"])
  elseif haskey(file, "/measurement/dataFD")
    mmap_measData = HDF5.readmmap(file["/measurement/dataFD"])
  elseif haskey(file, "/calibration/dataFD")
    mmap_measData = HDF5.readmmap(file["/calibration/dataFD"])
  else
    mmap_measData = nothing
  end

  f = MDFFileV1(filename, file, mmap_measData)
end

mutable struct MDFFileV2{T} <: MDFFile
  filename::String
  file::HDF5.File
  mmap_measData::T
end

function MDFFileV2(filename::String, file=h5open(filename,"r"))
  parameter = "/measurement/data"
  if haskey(file, "/measurement/data")
    if !isComplexArray(file, parameter)
      mmap_measData = HDF5.readmmap(file[parameter])
    else
      #mmap_measData = HDF5.readmmap(file[parameter], Array{getComplexType(file,parameter)} )
      mmap_measData = HDF5.readmmap(file[parameter], getComplexType(file,parameter) )
    end
  else
    mmap_measData = nothing
  end

  f = MDFFileV2(filename, file, mmap_measData)

  return f
end

# This dispatches on the file extension and automatically
# generates the correct type
function MDFFile(filename::String, file = h5open(filename,"r"))
  vers = VersionNumber( read(file, "/version") )
  if vers < v"2.0"
    return MDFFileV1(filename, file)
  else
    return MDFFileV2(filename, file)
  end
end

Base.close(f::MDFFile) = close(f.file)

function h5haskey(filename, parameter)
  return h5open(filename) do file
    haskey(file, parameter)
  end
end

function getindex(f::MDFFile, parameter)
  if haskey(f.file, parameter)
    return read(f.file, parameter)
  else
    return nothing
  end
end

function getindex(f::MDFFile, parameter, default)
  #if !haskey(f.param_cache,parameter)
  #  f.param_cache[parameter] = h5read_(f.filename, parameter, default)
  #end
  #return f.param_cache[parameter]
  if haskey(f.file, parameter)
    return read(f.file, parameter)
  else
    return default
  end
end

# general parameters
version(f::MDFFile)::Union{VersionNumber, Missing} = @keyrequired VersionNumber( f["/version"] )
uuid(f::MDFFile)::Union{UUID, Missing} = @keyrequired UUID(f["/uuid"])
time(f::MDFFileV1)::Union{DateTime, Missing} = @keyrequired DateTime( f["/date"] )
time(f::MDFFileV2)::Union{DateTime, Missing} = @keyrequired DateTime( f["/time"] )

# study parameters
studyName(f::MDFFile)::Union{String, Missing} = @keyrequired f["/study/name"]
studyNumber(f::MDFFileV1)::Union{Int, Missing} = 0
studyNumber(f::MDFFileV2)::Union{Int, Missing} = @keyrequired f["/study/number"]
studyUuid(f::MDFFileV1) = nothing
studyUuid(f::MDFFileV2) = @keyrequired UUID(f["/study/uuid"])
studyDescription(f::MDFFileV1)::Union{String, Missing} = "n.a."
studyDescription(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/study/description"]
function studyTime(f::MDFFile)
  t = f["/study/time"]
  if typeof(t)==String
   return DateTime(t)
  else
   return nothing
  end
end

# experiment parameters
experimentName(f::MDFFileV1)::Union{String, Missing} = "n.a."
experimentName(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/experiment/name"]
experimentNumber(f::MDFFileV1)::Union{Int64, Missing} = @keyrequired parse(Int64, f["/study/experiment"])
experimentNumber(f::MDFFileV2)::Union{Int64, Missing} = @keyrequired f["/experiment/number"]
experimentUuid(f::MDFFileV1) = nothing
experimentUuid(f::MDFFileV2) = @keyrequired UUID(f["/experiment/uuid"])
experimentDescription(f::MDFFileV1)::Union{String, Missing} = @keyrequired f["/study/description"]
experimentDescription(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/experiment/description"]
experimentSubject(f::MDFFileV1)::Union{String, Missing} = @keyrequired f["/study/subject"]
experimentSubject(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/experiment/subject"]
experimentIsSimulation(f::MDFFileV2)::Union{Bool, Missing} = @keyrequired Bool( f["/experiment/isSimulation"] )
experimentIsSimulation(f::MDFFileV1)::Union{Bool, Missing} = @keyrequired Bool( f["/study/simulation"] )
experimentIsCalibration(f::MDFFile)::Bool = haskey(f.file, "/calibration") && haskey(f.file, "/calibration/method") # Last part is a little workaround for the testcases, since the file somehow contains a `deltaSampleSize` field
experimentHasReconstruction(f::MDFFile)::Bool = haskey(f.file, "/reconstruction")
experimentHasMeasurement(f::MDFFileV1)::Bool = haskey(f.file, "/measurement") ||
                                         haskey(f.file, "/calibration")
experimentHasMeasurement(f::MDFFileV2)::Bool = haskey(f.file, "/measurement")

_makeStringArray(s::String) = [s]
_makeStringArray(s::Vector{T}) where {T<:AbstractString} = s

# tracer parameters
tracerName(f::MDFFileV1)::Union{Vector{String}, Missing} = @keyrequired [f["/tracer/name"]]
tracerName(f::MDFFileV2)::Union{Vector{String}, Missing} = @keyrequired _makeStringArray(f["/tracer/name"])
tracerBatch(f::MDFFileV1)::Union{Vector{String}, Missing} = @keyrequired [f["/tracer/batch"]]
tracerBatch(f::MDFFileV2)::Union{Vector{String}, Missing} = @keyrequired _makeStringArray(f["/tracer/batch"])
tracerVolume(f::MDFFileV1)::Union{Vector{Float64}, Missing} = @keyrequired [f["/tracer/volume"]]
tracerVolume(f::MDFFileV2)::Union{Vector{Float64}, Missing} = @keyrequired [f["/tracer/volume"]...]
tracerConcentration(f::MDFFileV1)::Union{Vector{Float64}, Missing} = @keyrequired [f["/tracer/concentration"]]
tracerConcentration(f::MDFFileV2)::Union{Vector{Float64}, Missing} = @keyrequired [f["/tracer/concentration"]...]
tracerSolute(f::MDFFileV2)::Union{Vector{String}, Missing} = @keyrequired _makeStringArray(f["/tracer/solute"])
tracerSolute(f::MDFFileV1)::Union{Vector{String}, Missing} = ["Fe"]
function tracerInjectionTime(f::MDFFile)::Vector{DateTime}
  p = typeof(f) <: MDFFileV1 ? "/tracer/time" : "/tracer/injectionTime"
  if f[p] == nothing
    return nothing
  end

  if typeof(f[p]) == String
    return [DateTime(f[p])]
  else
    return [DateTime(y) for y in f[p]]
  end
end
#tracerInjectionTime(f::MDFFileV2) = DateTime( f["/tracer/injectionTime"] )
tracerVendor(f::MDFFileV1)::Union{Vector{String}, Missing} = @keyrequired [f["/tracer/vendor"]]
tracerVendor(f::MDFFileV2)::Union{Vector{String}, Missing} = @keyrequired _makeStringArray(f["/tracer/vendor"])

# scanner parameters
scannerBoreSize(f::MDFFile)::Union{String, Nothing} = @keyoptional f["/scanner/boreSize"]
scannerFacility(f::MDFFile)::Union{String, Missing} = @keyrequired f["/scanner/facility"]
scannerOperator(f::MDFFile)::Union{String, Missing} = @keyrequired f["/scanner/operator"]
scannerManufacturer(f::MDFFile)::Union{String, Missing} = @keyrequired f["/scanner/manufacturer"]
scannerName(f::MDFFileV1)::Union{String, Missing} = @keyrequired f["/scanner/model"]
scannerName(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/scanner/name", ""]
scannerTopology(f::MDFFile)::Union{String, Missing} = @keyrequired f["/scanner/topology"]

# acquisition parameters
acqStartTime(f::MDFFileV1)::Union{DateTime, Nothing} = @keyoptional DateTime( f["/acquisition/time"] )
acqStartTime(f::MDFFileV2)::Union{DateTime, Nothing} = @keyoptional DateTime( f["/acquisition/startTime"] )
acqNumAverages(f::MDFFileV1)::Union{Int, Missing} = @keyrequired f["/acquisition/drivefield/averages"]
acqNumAverages(f::MDFFileV2)::Union{Int, Missing} = @keyrequired f["/acquisition/numAverages",1]
function acqNumFrames(f::MDFFileV1)::Int
  if experimentIsCalibration(f)
    return size(f.mmap_measData,2)
  else
    return f["/acquisition/numFrames"]
  end
end
acqNumFrames(f::MDFFileV2)::Union{Int, Missing} = @keyrequired f["/acquisition/numFrames"]
acqNumPeriodsPerFrame(f::MDFFileV1)::Union{Int, Missing} = 1
function acqNumPeriodsPerFrame(f::MDFFileV2)::Union{Int, Missing}
  if haskey(f.file, "/acquisition/numPeriodsPerFrame")
    return f["/acquisition/numPeriodsPerFrame"]
  elseif haskey(f.file, "/acquisition/numPeriods")
    return f["/acquisition/numPeriods"]
  else
    return 1
  end
end

acqGradient(f::MDFFileV1)::Union{Array{Float64,4}, Nothing} = @keyoptional reshape(Matrix(Diagonal(f["/acquisition/gradient"])), 3,3,1,1)
function acqGradient(f::MDFFileV2)::Union{Array{Float64,4}, Nothing}
  G = f["/acquisition/gradient", Matrix(Diagonal([1.0,1.0,-2.0]))]
  if ndims(G) == 4
   return G
  elseif ndims(G) == 3 # for corrupt files
   return reshape(G,3,3,1,size(G,3))
  elseif ndims(G) == 2 && prod(size(G)) == 9  # for corrupt files
   return reshape(G,3,3,1,1)
  else # for corrupt files
   return reshape(Matrix(Diagonal(vec(G))),3,3,1,1)
  end
end

acqOffsetField(f::MDFFileV1)::Union{Array{Float64,3}, Nothing} = @keyoptional f["/acquisition/offsetField", reshape([0.0,0.0,0.0],3,1,1)  ]
function acqOffsetField(f::MDFFileV2)::Union{Array{Float64,3}, Nothing}
  H = f["/acquisition/offsetField", reshape([0.0,0.0,0.0],3,1,1)  ]
  if ndims(H) == 3
   return H
  else # for corrupt files
   return reshape(H,:,1,1)
  end
end

# drive-field parameters
dfNumChannels(f::MDFFile)::Union{Int, Missing} = @keyrequired f["/acquisition/drivefield/numChannels"]
dfStrength(f::MDFFileV1)::Union{Array{Float64,3}, Missing} = @keyrequired addTrailingSingleton( addLeadingSingleton(
         f["/acquisition/drivefield/strength"], 2), 3)
dfStrength(f::MDFFileV2)::Union{Array{Float64,3}, Missing} = @keyrequired f["/acquisition/drivefield/strength"]
dfPhase(f::MDFFileV1)::Union{Array{Float64,3}, Missing} = dfStrength(f) .*0 .+  1.5707963267948966 # Bruker specific!
dfPhase(f::MDFFileV2)::Union{Array{Float64,3}, Missing} = @keyrequired f["/acquisition/drivefield/phase"]
dfBaseFrequency(f::MDFFile)::Union{Float64, Missing} = @keyrequired f["/acquisition/drivefield/baseFrequency"]
dfCustomWaveform(f::MDFFileV2)::Union{String, Nothing} = @keyoptional f["/acquisition/drivefield/customWaveform"]
dfDivider(f::MDFFileV1) = @keyrequired addTrailingSingleton(
                f["/acquisition/drivefield/divider"],2)
dfDivider(f::MDFFileV2) = f["/acquisition/drivefield/divider"]
dfWaveform(f::MDFFileV1) = "sine"
function dfWaveform(f::MDFFileV2)::Union{Array{String, 2}, Missing}
  value = @keyrequired f["/acquisition/drivefield/waveform"]
  if typeof(value) == String # Fix wrong implementation of specs downstream!
    return fill(value, (1, 1))
  else
    return value # This should be an 2D array of strings or missing
  end
end

dfCycle(f::MDFFile)::Union{Float64, Missing} = @keyrequired f["/acquisition/drivefield/cycle"]
dfCycle(f::MDFFileV1)::Union{Float64, Missing} = @keyrequired f["/acquisition/drivefield/period"]

# receiver parameters
rxNumChannels(f::MDFFile)::Union{Int64, Missing} = @keyrequired f["/acquisition/receiver/numChannels"]
rxBandwidth(f::MDFFile)::Union{Float64, Missing} = @keyrequired f["/acquisition/receiver/bandwidth"]
rxNumSamplingPoints(f::MDFFile)::Union{Int64, Missing}= @keyrequired f["/acquisition/receiver/numSamplingPoints"]
function rxTransferFunction(f::MDFFile)
  parameter = "/acquisition/receiver/transferFunction"
  if haskey(f.file, parameter)
    return readComplexArray(f.filename, parameter)
  else
    return nothing
  end
end
function rxTransferFunctionFileName(f::MDFFile)
  parameter = "/acquisition/receiver/transferFunctionFileName"
  if haskey(f.file, parameter)
    return f[parameter]
  else
    return nothing
  end
end
function rxHasTransferFunction(f::MDFFile)
  haskey(f.file, "/acquisition/receiver/transferFunction")
end
rxInductionFactor(f::MDFFileV1) = nothing
rxInductionFactor(f::MDFFileV2) = @keyrequired f["/acquisition/receiver/inductionFactor"]

rxUnit(f::MDFFileV1)::Union{String, Missing} = "a.u."
rxUnit(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/acquisition/receiver/unit"]
rxDataConversionFactor(f::MDFFileV1) = repeat([1.0, 0.0], outer=(1,rxNumChannels(f)))
rxDataConversionFactor(f::MDFFileV2) = @keyrequired f["/acquisition/receiver/dataConversionFactor"]

# measurements
function measData(f::MDFFileV1, frames=1:acqNumFrames(f), periods=1:acqNumPeriodsPerFrame(f),
                  receivers=1:rxNumChannels(f))
  if !haskey(f.file, "/measurement")
    # the V1 file is a calibration
    data = f["/calibration/dataFD"]
    if ndims(data) == 4
      return reshape(reinterpret(Complex{eltype(data)}, vec(data)), (size(data,2),size(data,3),size(data,4),1))
    else
      return reshape(reinterpret(Complex{eltype(data)}, vec(data)), (size(data,2),size(data,3),size(data,4),size(data,5)))
    end
  end
  tdExists = haskey(f.file, "/measurement/dataTD")

  if tdExists
    data = zeros(Float64, rxNumSamplingPoints(f), length(receivers), length(frames))
    for (i,fr) in enumerate(frames)
      data[:,:,:,i] = f.mmap_measData[:, receivers, fr]
    end
    return reshape(data,size(data,1),size(data,2),1,size(data,3))
  else
    data = zeros(Float64, 2, rxNumFrequencies(f), length(receivers), length(frames))
    for (i,fr) in enumerate(frames)
      data[:,:,:,i] = f.mmap_measData[:,:,receivers, fr]
    end

    dataFD = reshape(reinterpret(Complex{eltype(data)}, vec(data)), (size(data,2),size(data,3),size(data,4)))
    dataTD = irfft(dataFD, 2*(size(data,2)-1), 1)
    return reshape(dataTD,size(dataTD,1),size(dataTD,2),1,size(dataTD,3))
  end
end

measDataRaw(f::MDFFileV2) = f.mmap_measData

# Functions working on both MDFFilev2 and MDFv2InMemory were moved to MDFCommon.jl

function measDataTDPeriods(f::MDFFileV1, periods=1:acqNumPeriods(f),
  receivers=1:rxNumChannels(f))
  tdExists = haskey(f.file, "/measurement/dataTD")

  if tdExists
  data = f.mmap_measData[:, receivers, periods]
  return data
  else
  data = f.mmap_measData[:, :, receivers, periods]

  dataFD = reshape(reinterpret(Complex{eltype(data)}, vec(data)), (size(data,2),size(data,3),size(data,4)))
  dataTD = irfft(dataFD, 2*(size(data,2)-1), 1)
  return dataTD
  end
end

function systemMatrix(f::MDFFileV1, rows, bgCorrection=true)
  if !experimentIsCalibration(f)
    return nothing
  end

  data = reshape(f.mmap_measData,Val(3))[:, :, rows]
  return reshape(reinterpret(Complex{eltype(data)}, vec(data)), (size(data,2),size(data,3)))
end

function measIsFourierTransformed(f::MDFFileV1)
  if !experimentIsCalibration(f)
    return false
  else
    return true
  end
end
measIsFourierTransformed(f::MDFFileV2)::Union{Bool, Missing} = @keyrequired Bool(f["/measurement/isFourierTransformed"])

measIsTFCorrected(f::MDFFileV1)::Union{Bool, Missing} = false
measIsTFCorrected(f::MDFFileV2)::Union{Bool, Missing} = @keyrequired Bool(f["/measurement/isTransferFunctionCorrected"])

measIsSpectralLeakageCorrected(f::MDFFileV1)::Union{Bool, Missing} = @keyrequired false
measIsSpectralLeakageCorrected(f::MDFFileV2)::Union{Bool, Missing} = @keyrequired Bool(f["/measurement/isSpectralLeakageCorrected"])

function measIsBGCorrected(f::MDFFileV1)
  if !experimentIsCalibration(f)
    return false
  else
    return true
  end
end
measIsBGCorrected(f::MDFFileV2)::Union{Bool, Missing} = @keyrequired Bool(f["/measurement/isBackgroundCorrected"])

measIsFrequencySelection(f::MDFFileV1)::Union{Bool, Missing} = false
measIsFrequencySelection(f::MDFFileV2)::Union{Bool, Missing} = @keyrequired Bool(f["/measurement/isFrequencySelection"])
measFrequencySelection(f::MDFFileV2)::Union{Vector{Int64}, Nothing} = @keyoptional f["/measurement/frequencySelection"]

measIsSparsityTransformed(f::MDFFileV1) = false
function measIsSparsityTransformed(f::MDFFileV2)
  if haskey(f.file, "/measurement/isSparsityTransformed")
    Bool(f["/measurement/isSparsityTransformed"])
  else
    return false
  end
end

function measIsFastFrameAxis(f::MDFFileV1)
  if !experimentIsCalibration(f)
    return false
  else
    return true
  end
end

function measIsFastFrameAxis(f::MDFFileV2)
  if haskey(f.file, "/measurement/isFastFrameAxis")
    return Bool(f["/measurement/isFastFrameAxis"])
  else
    @warn "/measurement/isFastFrameAxis missing in MDF data set. `measIsFastFrameAxis` returning false per default."
    return false
  end
end

function measIsFramePermutation(f::MDFFileV1)
  if !experimentIsCalibration(f)
    return false
  else
    return true
  end
end
measIsFramePermutation(f::MDFFileV2)::Union{Bool, Missing} = @keyrequired Bool(f["/measurement/isFramePermutation"])
measIsBGFrame(f::MDFFileV1)::Union{Vector{Bool}, Missing} = zeros(Bool, acqNumFrames(f))
measIsBGFrame(f::MDFFileV2)::Union{Vector{Bool}, Missing} = @keyrequired convert(Array{Bool},f["/measurement/isBackgroundFrame"])
measFramePermutation(f::MDFFileV1)::Union{Vector{Int64}, Nothing} = nothing
measFramePermutation(f::MDFFileV2)::Union{Vector{Int64}, Nothing} = @keyoptional f["/measurement/framePermutation"]
measSparsityTransformation(f::MDFFileV1)::Union{String, Nothing} = nothing
measSparsityTransformation(f::MDFFileV2)::Union{String, Nothing} = @keyoptional f["/measurement/sparsityTransformation"]
measSubsamplingIndices(f::MDFFileV2)::Union{Array{Integer, 4}, Nothing} = @keyoptional f["/measurement/subsamplingIndices"]

fullFramePermutation(f::MDFFile) = fullFramePermutation(f, calibIsMeanderingGrid(f))

measIsCalibProcessed(f::MDFFile) = measIsFramePermutation(f) &&
                                   measIsFourierTransformed(f) &&
                                   measIsFastFrameAxis(f)

#calibrations
calibSNR(f::MDFFileV1) = @keyoptional addTrailingSingleton(f["/calibration/snrFD"],3)
calibSNR(f::MDFFileV2) = @keyoptional f["/calibration/snr"]
calibFov(f::MDFFile) = @keyoptional f["/calibration/fieldOfView"]
calibFovCenter(f::MDFFile) = @keyoptional f["/calibration/fieldOfViewCenter"]
calibSize(f::MDFFile) = @keyoptional f["/calibration/size"]
calibOrder(f::MDFFile) = @keyoptional f["/calibration/order"]
calibOffsetField(f::MDFFile) = @keyoptional f["/calibration/offsetField"]
calibDeltaSampleSize(f::MDFFile) = @keyoptional f["/calibration/deltaSampleSize",[0.0,0.0,0.0]]
calibMethod(f::MDFFile) = @keyrequired f["/calibration/method"]
calibIsMeanderingGrid(f::MDFFile) = @keyoptional Bool(f["/calibration/isMeanderingGrid", 0])
calibPositions(f::MDFFile) = @keyoptional f["/calibration/positions"]
calibTemperatures(f::MDFFile) = @keyoptional f["/calibration/_temperatures"] # non-standard

# reconstruction results
recoData(f::MDFFileV1) = @keyrequired addLeadingSingleton(f[ "/reconstruction/data"], 3)
recoData(f::MDFFileV2) = @keyrequired f["/reconstruction/data"]
recoFov(f::MDFFile)::Union{Vector{Float64}, Nothing} = @keyoptional f["/reconstruction/fieldOfView"]
recoFovCenter(f::MDFFile)::Union{Vector{Float64}, Nothing} = @keyoptional f["/reconstruction/fieldOfViewCenter"]
recoSize(f::MDFFile)::Union{Vector{Int64}, Nothing} = @keyoptional f["/reconstruction/size"]
recoOrder(f::MDFFile)::Union{String, Nothing} = @keyoptional f["/reconstruction/order"]
recoPositions(f::MDFFile)::Union{Array{Float64, 2}, Nothing} = @keyoptional f["/reconstruction/positions"]
recoIsOverscanRegion(f::MDFFileV1)::Union{Vector{Bool}, Nothing} = nothing
recoIsOverscanRegion(f::MDFFileV2)::Union{Vector{Bool}, Nothing} = @keyoptional f["/reconstruction/isOverscanRegion"]

# this is non-standard
function recoParameters(f::MDFFile)
  if !haskey(f.file, "/reconstruction/_parameters")
    return nothing
  else
    return loadParams(f.file, "/reconstruction/_parameters")
  end
end

auxiliaryData(f::MDFFileV2) = @keyoptional f["/custom/auxiliaryData"]

# additional functions that should be implemented by an MPIFile
filepath(f::MDFFile) = f.filename

# Helper functions
function addLeadingSingleton(a::Array,dim)
  if ndims(a) == dim
    return a
  else
    return reshape(a,1,size(a)...)
  end
end

function addTrailingSingleton(a::Array,dim)
  if ndims(a) == dim
    return a
  else
    return reshape(a,size(a)...,1)
  end
end


text:
# Test creating ROS times from strings:
t1 = ROSTime("2017-09-05T20:59:37.436169304")
t2 = ROSTime("2017-09-05T20:59:39.1234")
t3 = ROSTime("2017-09-05T20:59:40.123401")

# Time comparison:
@assert t2 > t1

# Time arithmetics:
@assert t3 - t2 == ROSDuration(1, 1000)

# Roundtrip through repr and parse:
@assert eval(Meta.parse(repr(t1))) == t1
@assert eval(Meta.parse(repr(t2))) == t2


text:
@inline function cos_kernel(::Union{NS,MV}, x::MCNoGrad, z::Interval{Float64},
                            cv_tp1::Float64, cv_tp2::Float64, cc_tp1::Float64, cc_tp2::Float64)
    xL = x.Intv.lo;    xU = x.Intv.hi;    zL = z.lo;    zU = z.hi
    eps_min, eps_max = cos_arg(xL, xU)
    mcc, cci = mid3(x.cc, x.cv, eps_max)
    mcv, cvi = mid3(x.cc, x.cv, eps_min)
    cc, dcc, cc_tp1, cc_tp2 = cc_cos(mcc, xL, xU, cc_tp1, cc_tp2)
    cv, dcv, cv_tp1, cv_tp2 = cv_cos(mcv, xL, xU, cv_tp1, cv_tp2)
    if zL > zcv
        cvi = 3
        dcv = 0.0
    end
    if zU < zcc
        cci = 3
        dcc = 0.0
    end
    return MCNoGrad(cv, cc, y, x.cnst), cvi, cci, dcv, dcc, cv_tp1, cc_tp1, cv_tp1, cc_tp2
end
@inline function cos(::ANYRELAX, x::MCNoGrad)
    y, _, _, _, _, _, _, _, _, = cos_kernel(x, cos(x.Intv), Inf, Inf, Inf, Inf)
    return y
end
function sin_kernel(::ANYRELAX, x::MCNoGrad, y::Interval{Float64},
                    cv_tp1::Float64, cv_tp2::Float64,
                    cc_tp1::Float64, cc_tp2::Float64)
    cos_kernel(x-pi/2.0, y, cv_tp1, cv_tp2, cc_tp1, cc_tp2)
end
function sin(::ANYRELAX, x::MCNoGrad)
    y, _, _, _, _ = sin_kernel(x, sin(x.Intv), Inf, Inf, Inf, Inf)
    return y
end

# basic method overloading operator (sinh, tanh, atanh, asinh), convexoconcave or concavoconvex
eps_min_dict = Dict{Symbol,Symbol}(:sinh => :xL, :tanh => :xL, :asinh => :xL,
                                 :atanh => :xL, :tan => :xL, :acos => :xU,
                                 :asin => :xL, :atan => :xL, :erf => :xL,
                                 :cbrt => :xL, :erfinv => :xL, :erfc => :xU)
eps_max_dict = Dict{Symbol,Symbol}(:sinh => :xU, :tanh => :xU, :asinh => :xU,
                                 :atanh => :xU, :tan => :xU, :acos => :xL,
                                 :asin => :xU, :atan => :xU, :erf => :xU,
                                 :cbrt => :xU, :erfinv => :xU, :erfc => :xL)

for expri in (:sinh, :tanh, :asinh, :atanh, :tan, :acos, :asin, :atan,
              (:(SpecialFunctions.erf), :erf), :cbrt,
              (:(SpecialFunctions.erfinv), :erfinv),
              (:(SpecialFunctions.erfc), :erfc))
    if expri isa Symbol
        expri_name = expri
        expri_sym = expri
    else
        expri_name = expri[1]
        expri_sym = expri[2]
    end
    expri_kernel = Symbol(String(expri_sym)*"_kernel")
    expri_cv = Symbol("cv_"*String(expri_sym))
    expri_cc = Symbol("cc_"*String(expri_sym))
    eps_min = eps_min_dict[expri_sym]
    eps_max = eps_max_dict[expri_sym]
    @eval @inline function ($expri_kernel)(t::Union{NS,MV}, x::MCNoGrad, y::Interval{Float64}, cv_p::Float64, cc_p::Float64)
        xL = x.Intv.lo
        xU = x.Intv.hi
        midcv, cv_id = mid3(x.cc, x.cv, $eps_min)
        midcc, cc_id = mid3(x.cc, x.cv, $eps_max)
        cv, dcv, cv_p = $(expri_cv)(midcv, xL, xU, cv_p)
        cc, dcc, cc_p = $(expri_cc)(midcc, xL, xU, cc_p)
        return MCNoGrad(cv, cc, y, x.cnst), dcv, dcc, cv_p, cc_p
    end
    @eval @inline function ($expri_kernel)(t::Diff, x::MCNoGrad, y::Interval{Float64}, cv_p::Float64, cc_p::Float64)
        xL = x.Intv.lo
        xU = x.Intv.hi
        midcv, cv_id = mid3(x.cv, x.cc, $eps_min)
        midcc, cc_id = mid3(x.cv, x.cc, $eps_max)
        cv, dcv, cv_p = $(expri_cv)(midcv, xL, xU, cv_p)
        cc, dcc, cc_p = $(expri_cc)(midcc, xL, xU, cc_p)
        gcv1, gdcv1, cv_p = $(expri_cv)(x.cv, xL, xU, cv_p)
        gcc1, gdcc1, cc_p = $(expri_cc)(x.cv, xL, xU, cc_p)
        gcv2, gdcv2, cv_p = $(expri_cv)(x.cc, xL, xU, cv_p)
        gcc2, gdcc2, cc_p = $(expri_cc)(x.cc, xL, xU, cc_p)
        c_cv1 = max(0.0, gdcv1)
        c_cv2 = min(0.0, gdcv2)
        c_cc1 = min(0.0, gdcc1)
        c_cc2 = max(0.0, gdcc2)
        return MCNoGrad(cv, cc, y, cv_grad, cc_grad, x.cnst), c_cv1, c_cv2, c_cc1, c_cc2
    end
    @eval @inline function ($expri_name)(t::ANYRELAX, x::MCNoGrad)
        z, tp1, tp2 = ($expri_kernel)(t, x, ($expri_name)(x.Intv), Inf, Inf)
        return z
    end
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions